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INTRODUCTION  TO  DIMSIMS 


LITERATURE  SURVEY 

Diagonally  implicit  multistage  integration  methods  (DIMSIMs)  were  first  described  in  a 
1992  paper  by  John  Butcher  (Reference  1)  in  which  he  laid  out  the  essential  elements  of  a 
new  family  of  general  linear  methods.  His  stated  purpose  was  to  overcome  the  glaring 
weaknesses  of  existing  methods,  that  is,  lack  of  A-stability  for  high  order  linear  multistep 
methods  and  low  stage  order  and  high  implementation  costs  for  A-stable  implicit  Runge- 
Kutta  methods.  These  methods  are  diagonally  implicit  and  hence  have  computational 
complexity  properties  similar  to  diagonally  implicit  Runge-Kutta  (DIRK)  methods,  but 
utilize  additional  parameters  generated  through  the  general  linear  design  to  overcome  the 
stage  order  (and  hence  stiff  order)  limitations  of  DIRK  methods.  Explicit  DIMSIMs  are  not 
technically  “diagonally  implicit”,  since  the  diagonal  is  actually  0,  but  have  an  advantage 
over  Runge-Kutta  methods  in  that  the  order  barriers  for  explicit  Runge-Kutta  methods  do 
not  apply,  and  p-stage  methods  of  order  p  do  indeed  exist  for  all  positive  integers  p.  The 
original  concept  called  for  stage  order  q  to  equal  the  number  of  internal  stages  s,  the 
number  of  external  stages  r,  and  the  order  p.  In  a  subsequent  paper  with  Jackiewicz 
(Reference  2)  the  family  of  DIMSIMs  was  extended  to  include  adjacent  methods  for  which 
s+l=r  =  q,  p  =  qorq+l,s  =  r+  l  =  q,  p=qorq+l,  and  s  =  r  =  q,  p  =  q+l.  A 
significant  step  toward  practical  utilization  was  taken  with  another  paper  by  Butcher  and 
Jackiewicz  (Reference  3)  that  lays  out  techniques  for  error  estimation,  interpolation,  and 
step-size  changing.  Jackiewicz,  Vermiglio,  and  Zennaro  (Reference  4)  devised  an 
alternative  step-size  changing  strategy  and  showed  how  incorporation  of  an  additional 
external  stage  could  provide  a  satisfactory  continuous  method.  In  a  separate  paper 
(Reference  5)  they  also  showed  that  there  exist  explicit  DIMSIMs  with  regularity  properties 
not  possessed  by  explicit  Runge-Kutta  methods.  Butcher,  Chartier,  and  Jackiewicz,  in  an 
unpublished  manuscript,  “Nordsieck  representation  of  DIMSIMs,”  recently  proposed  an 
alternative  representation  of  DIMSIMs  with  promise  of  simplifying  analysis  and 
implementation.  Both  explicit  and  implicit  DIMSIMs  up  to  the  order  8  have  now  been 
found  with  appropriate  stability  properties  and  were  announced  by  Butcher,  Jackiewicz, 
and  Mittelmann  (Reference  6),  extending  techniques  described  in  earlier  papers  by  Butcher 
and  Jackiewicz  (Reference  7  and  8). 


BASIC  DEFINITIONS  AND  RELATIONSHIPS 


We  consider  the  initial  value  problem: 
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We  define  matrices  A,  U,  B  and  V  such  that  A  is  s  x  s,  V  is  r  x  r,  U  is  s  x  r,  and  B  is  r  x  s . 
Let  Y  be  composed  of  the  s  internal  stages,  F  be  composed  of  the  s  stage  derivatives 

( Fj  =  f{tn  +  Cjh,  F7  j)  and  be  composed  of  the  r  external  stage  values.  Then  if  h  is  the 
step  size,  the  solution  advances  one  step  through  the  relationships: 

Yi  =  h  I  OijFj  +  X  Uijy^~l\  i  =  1,2,..., s 
7-1  7=1 

y\n]  =hi  bijFj  +  £  Vy-yM,  i  =  1,2, ...,r  (2) 

7=1  7=1 

n-\,2,...,N,  Nh  =  T-t0 


The  external  stages  are  defined  through  Taylor  expansion  so  that  if 

7=0 


(3) 


then  we  must  have,  for  some  constants  , 

y{"'=  ia9Vyw(»„)  +  0(A^1), 

7=0 


(4) 


for  a  method  of  order  p.  Butcher  has  observed  (Reference  9)  that  the  effect  then  is  to 
calculate  r  neighboring  trajectories  and  to  use  these  to  determine  solution  values.  The  s 
values  C;  are  chosen  initially  and  other  method  parameters  are  then  determined  in  a  way  to 
produce  high  stage  order,  that  is,  so  that: 


Yt  =  y{tn_i  +  hci )  +  0{hq+l),  i  =  1,2,. ..,5. 


(5) 


where  q  is  defined  to  be  the  stage  order. 

The  description  “diagonally  implicit”  comes  from  the  form  of  the  matrix  A,  and  the 
following  discussion  is  true  for  both  DIMSIMs  and  DIRK  methods.  If  we  restrict  A  to  be 
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lower  triangular  with  a  constant  along  the  diagonal,  the  complexity  of  the  solution  of  the 
system  of  nonlinear  equations  determining  Y  in  Equation  2  reduces  greatly.  If  A  is  dense 
and  a  standard  gaussian  elimination  approach  is  used  in  a  modified  Newton  method,  the 
arithmetic  complexity  is  0((Ms)3),  or  0(M3s3),  where  s  is  the  number  of  stages  of  the 
method  and  M  is  the  dimension  of  the  initial  value  problem.  Simply  requiring  A  to  be 
lower  triangular  separates  the  system  of  Ms  simultaneous  nonlinear  equations  into  s 
systems  of  M  simultaneous  equations  to  be  solved  in  sequence,  each  with  arithmetic 
complexity  0(M3)  for  total  complexity  0(sM3),  a  reduction  by  a  factor  of  0(s2). 

The  advantage  of  having  a  single  value  along  the  diagonal  may  be  seen  from  a  closer 
examination  of  the  nonlinear  system  solution  process,  which  typically  involves  a  modified 
Newton  method.  The  Newton  iteration  to  solve  the  nonlinear  system  g(x)  =  0  takes  the 
form 


xn+ 1  ~  xn  J{xn)  §(Xn)’ 


where  J  is  the  Jacobian  of  g.  Of  course  the  inverse  of  the  Jacobian  is  not  actually  calculated 
and  instead  a  technique  such  as  gaussian  elimination  is  utilized,  and  the  process  followed  is 
to  calculate  a  correction: 


Jr(xrt)^n+ 1  s{xn)’ 

xn+\  =  xn  +  <W 


An  LU  or  PLU  factorization  of  J(xn)  is  called  for  here  at  each  iteration,  which  is  0(M3), 
where  M  is  the  size  of  the  system,  and  this  is  the  most  time  consuming  step.  In  practice  a 
new  Jacobian  is  evaluated  and  a  new  matrix  factorization  is  carried  out  only  when 
convergence  seems  too  slow.  In  the  case  of  DIMSIMs  the  equation  for  Y  takes  the  form 


Y1=hTaJtf{tn.l+cth,rk)  +  hajif(l„^  +Cjh,Yj)  +  #-1\  (8) 


with  U=I  as  is  usually  the  case.  Then  we  are  solving  an  equation  of  the  form  gj(Yj)=0, 
where  g-  takes  the  form: 


gj(Yj)=rj -ItaA'n- 1  1  +ckh,Yt)-ylf-'\ 


(9) 


5 


NAWCWPNS  TP  8340 


All  the  dependency  on  Yj  is  contained  in  the  first  two  terms  and  so  only  these  terms  affect 
the  calculation  of  the  Jacobian.  Now  if  all  diagonal  elements  a-  are  the  same  the  Jacobians 
J  will  vary  from  stage  to  stage  only  with  the  solution,  and  new  Jacobians  and  their  LU 
decompositions  will  only  have  to  be  computed  when  convergence  seems  too  slow,  which 
should  be  relatively  rare  for  small  step  size  h.  This  can  result  in  substantial  savings. 

A  similar  reduction  follows  from  utilization  of  a  nonsingular  A  (perhaps  dense)  with  a 
single  eigenvalue,  and  this  has  been  applied  in  the  implicit  Runge-Kutta  code  STRIDE 
developed  by  Burrage,  Butcher,  and  Chipman  (Reference  10).  Because  of  the  work 
required  for  linear  transformation,  this  is  to  be  avoided  if  possible  due  to  the  number  of 
transformations  that  become  necessary.  Of  course  the  low  stage  order  of  DIRK  methods, 
which  drastically  reduces  the  observed  order  of  the  method  for  stiff  equations  to  second 
order,  made  this  alternative  approach  attractive  for  use  in  STRIDE. 

The  conditions  (Equations  2  through  5)  may  be  re-expressed  in  a  convenient  form  as 
follows.  Let  W  be  the  r  x  (p  +  1)  matrix  of  oc^  values,  and  we  denote  the  vector  consisting 
of  the  kth  column  of  W  as  a^.Let  Z  be  a  vector  with  element  Zj  =  zs'\  Then  define  w(z)  = 
WZ.  Furthermore  we  may  define  ecz  as  the  vector 


„CZ 


(10) 


Then  the  following  theorem  may  be  used  in  determining  the  coefficients  of  the  method: 
Theorem  1:  A  DIMSIM  (Equation  2)  has  order  p  and  stage  order  p  if  and  only  if 


ecz  =  zAecz  +  Uw(z)  +  o(zp+1^j, 
ezw(z )  =  zBecz  +  Vw(z)  +  o(zp+1). 


(11) 


Proof:  The  proof  is  given  by  Butcher  (Reference  1)  and  is  included  here  for  the  sake  of 
completeness. 

The  stage  order  condition  gives 
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We  apply  the  derivative  function  f  to  both  sides  and  find  that 

/(<„-,  +  CM)  =  /(<„-!  +  +  cih)  +  O^1)) 

=  /('„-!  +cih,y(t„_1  +  c,/i))  +  o(/i'"1) 

We  then  may  write,  using  Taylor  expansion, 

¥(<„-i +c{h^). 

Since 

j= o 

we  find  through  Taylor  expansion  that 

y!"1  =  f  {  £  idy-iWl)*' 

;=ou=o*!  / 

We  now  utilize  Equation  2,  provide  a  Taylor  expansion  for  each  Yj,  and  use  W  to  express 
eachyt”-1]  in  terms  of  hky^k\tn_ j).  Then  we  can  combine  coefficients  to  obtain  the  two 
equations 
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p 

I 

k=0 


p 

1 

k=0 


4  -  ikaijC)-1  -  iuijaijkl\^4k\tn_1)  =  o(hp+1) 


j= 1  J 


( 


X  -  X  *Vf '  -  X  W!  =  °Kl) 

=0  l-  7=1  7=1  J  K- 


V 


Now  each  coefficient  (call  them  dk  and  dk,  respectively  for  the  two  equations)  of 


hk 


k\ 
have 


y^(?n_l)  (for  each  k)  may  be  set  to  0  since  they  are  all  independent  of  h.  Then  we 


p  z 
Z.dkZ-  =  0, 
k= 0  k" 

P  ~  zk 

14  7T  =  0^ 

k=0  «! 


which  leads  to  two  new  equations.  The  first  is 


p  .  p  s  ,  .  p  r  . 

X  C?J7-  X  Zto,j4 -'jr-  X  X«,^=  0. 

it=o  k!  k-o  7=i  k.  ^=o  7=i 


The  first  summation  may  be  identified  as  the  first  p  +  1  terms  in  the  Taylor  expansion  of 
ec,z,  leaving  a  difference  that  is  of  order  0(zp+1).  Reversing  the  order  of  the  summations, 

r 

the  third  term  may  be  readily  identified  as  -  'ZuijWj(z).  Reversing  the  order  of  summation 

7=1 

for  the  second  term  and  factoring  as  appropriate,  we  have 


J  p  ir-i  zk  1  s  p~l  k  zk 
I  ayZ  I  c  77  777  =  I  aijz  I  cj  T7- 

7= l  1  k= l  y  (k  - 1)!  7=1  *= o  k\ 
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c  ■ z 

Here  the  sum  over  k  may  be  identified  as  the  first  p  terms  in  the  Taylor  expansion  of  e  1  , 

leaving  a  difference  that  is  0(zp).  Then  the  second  term  is  -  'ZaijzeCjZ  +  o(zp+l ).  Thus 

7=1 

the  first  equation  is  equivalent  to: 


ec‘z  -  JtzageejZ  -  'LugWjiz)  =  o(zp+l), 
7=1  7=1 

or,  in  matrix  form, 


ecz  -  zAecz  +  Uw(z)  +  o[zp+^ . 


The  second  new  equation  is: 


p  k  b\  7k  p  s  ,  .  7*  p  r  7k 

I  X  —cCik-i- —  X  X*V/'  —  ■ -  x  lvijajkk\Zr  =  0. 

£o,to  /!  1  k\  £0jt i  tJ  J  k\  *tb mij  Jk  k\ 


We  interchange  the  order  of  summation  and  find  for  the  second  term  similarly  as  for  the 

second  term  above  that  we  have  -  ’EbijZeCjZ  +  o(z-p+1j.  For  the  third  term  our  expression 

7=1 

r 

is  similar  to  that  for  the  third  term  in  the  first  equation  and  we  have  -Xvl7w.(z) .  For  the 

,  7=1 

first  term  we  may  observe  that  it  is  e2wi(z)+0(zp+1)  by  looking  at  ebv^z)  as  the  product  of 
the  first  p+1  terms  in  the  Taylor  expansion  for  ez  times  the  terms  in  W;(z).  If  terms  of  the 
same  order  in  z  are  combined  and  terms  of  order  in  z  higher  than  p  are  dropped  the 
expressions  may  be  seen  to  be  identical.  Then  the  second  new  equation  may  be  rewritten 
as 


ezWi{z)~  Ibyze CjZ  -  'Lvijwj(z )  =  o(zp+l\ 

7=1  7=1 

or  in  matrix  form, 


ezw(z)  -  zBecz  -  Vecz  =  o[zp+l ) . 
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These  are  equivalent  to  the  second  equation  in  the  theorem.  On  the  other  hand,  these  steps 
may  be  reversed  to  obtain  the  order  and  stage  order  conditions.  ■ 

It  turns  out  that  a  restriction  on  the  coefficient  matrix  B  provides  assurance  that  order 
and  stage  order  conditions  are  met,  according  to  the  following  theorem  derived  by  Butcher 
(Reference  1).  The  symbol  e  will  be  used  throughout  this  report  for  the  vector  of  the 
appropriate  length  for  the  context  for  which  each  element  is  1. 

Theorem  2.  Let  r  =  s  =  p,  Ve  =  e.  Then  the  DIMSIM 


A  I 
B  V 


is  of  order  p  and  stage  order  q  =  p  if  and  only  if 


B=B0-AB,-VB2+VA, 


where 


(Pj{^  +  ci) 
^j{cj)  ’ 


&2  ,i,j  - 


_  \o<Pj(x)dx 

tj(cj) 


(12) 


and  where 


<Pj(x)=  Il(x-ck). 
k±j 


(13) 


Proof:  See  Butcher  (Reference  1). 

Note  that  using  this  theorem  eliminates  the  elements  of  B  as  free  parameters  when 
deriving  methods.  It  also  has  the  following  immediate  corollary,  since  for  any  specified 
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vector  c  and  matrices  A  and  V  a  construction  may  be  completed,  leaving  stability  as  the 
principal  issue  in  deriving  new  methods. 

Corollary:  For  each  integer  p  >  1,  DIMSIMs  of  order  p  and  stage  order  q  =  p  exist  for 
s  =  r  =  p,  U  =  I,  where  s  is  the  number  of  internal  stages  and  r  is  the  length  of  the  external 
stage  vector. 

In  what  follows  we  will  restrict  ourselves  to  DIMSIMs  for  which  s  =  r  =  p  =  q  and 
U=  I. 


STABILITY  AND  CONSISTENCY  PROPERTIES  OF  DIMSIMS 

We  may  note  that  in  order  to  handle  the  simple  scalar  equation  y'  =  0  we  must  have  the 
preconsistency  condition  Ve  =  e.  The  eigenvalues  of  V  detemine  the  power-boundedness 
of  the  method,  required  for  zero-stability,  and  these  must  all  have  magnitudes  not  greater 
than  1,  and  those  with  magnitude  equal  to  one  must  have  one-dimensional  Jordan  blocks. 
A  typical  design  choice  is  to  choose  all  eigenvalues  to  be  0  except  for  the  unit  eigenvalue 
associated  with  eigenvector  e  as  is  the  case  for  Runge-Kutta  methods,  leaving  V  of  rank  1 
and  all  rows  identical. 

The  standard  consistency  condition,  related  to  the  solution  of  the  equation  y'  =  1,  then 
requires  that  there  must  exist  a  consistency  vector  u  such  that  Be  +  Vu  =  e  +  u  (Reference 
1 1).  We  first  note  that  equating  the  terms  of  order  zero  in  z  in  the  Taylor  expansion  about 
z  =  0  for  the  first  equation  of  Equation  1 1  tells  us  that  a0j  =  1,  j  =  l,...,p,  so  that  a0  =  e 
(a0  is  the  first  column  of  W).  It  follows  that  Ve  =  e  and  Be  +  Va,  =  e  +  a,  where  a,  is  the 
second  column  of  W,  made  up  of  elements  alj;  as  may  be  seen  by  equating  zeroth  and  first 
order  terms  in  z  (respectively)  in  the  Taylor  expansion  about  z  =  0  for  the  second  equation 
of  Equation  1 1  and  hence  a,  is  a  consistency  vector  for  the  method. 

The  stability  matrix  for  a  DIMSIM  is  (Reference  1) 


M(z)  =  V  +  zB(I-zA)'\ 


(14) 


This  is  easily  seen  from  the  method  (Equation  2)  using  the  standard  test  problem  y’  =  A,y, 
y(t0)  =  y0.  We  solve  the  first  equation  in  2  for  Y  and,  noting  that  F(Y)  =  XY,  we  obtain, 
setting  z  =  hX, 


y=(/-zArV"-11. 


Then  we  may  substitute  this  expression  in  the  second  equation  of  2  to  obtain 
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y[n]  =  (zB{I  -  zA)~l  +  y)j[n"i:i. 


Thus  the  region  of  absolute  stability  of  a  DIMSIM  is  the  region 
S  =  [z :  we  cr(M(z))  =>  w  <  lj. 

If  S  includes  the  entire  open  left  half  plane  the  method  is  called  A-stable.  We  also  define 
the  associated  stability  polynomial 


p(co,z)  =  det(w/  -  M(z)) . 


A  method  may  typically  be  verified  to  be  A-stable  either  by  using  the  Schur  criterion  (see, 
for  example,  Lambert  (Reference  13)),  or  by  reducing  the  stability  polynomial  to  a  familiar 
form  associated  with  a  Runge-Kutta  method  known  to  be  A-stable. 


A  FEW  SIMPLE  EXAMPLES  OF  DIMSIMS 


Butcher  divides  DIMSIMs  into  4  categories,  or  types,  using  criteria  of  explicit  or 
implicit  and  suitable  or  not  suitable  for  parallel  evaluation  of  stages.  Note  that  for  a  matrix 
A  that  is  lower  triangular  as  is  required  for  DIMSIMs,  if  it  is  also  diagonal  (including  the 
zero  matrix)  the  stage  evaluations  are  completely  independent  of  one  another  and  may  be 
carried  out  in  parallel.  Thus  we  have  the  following  order  2  examples  with  c  =  [0,1],  first 
developed  by  Butcher  and  illustrating  his  taxonomy  (Reference  1): 


Type  1  (explicit,  serial): 


'A  U~ 
B  V 


'00  1  O' 

2  0  0  1 

1111 
4  4  2  2 

2  _1  1  1 
4  4  2  2 


(15) 


This  method  has  by  design  the  same  stability  function  (and  region)  as  a  2-stage  explicit 
Runge-Kutta  methods  of  order  2.  That  is,  the  eigenvalues  of  M(z)  are  zero  and  R(z)  where 
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the  stability  function  of  a  Runge-Kutta  method  of  order  2. 
Type  2  (implicit,  serial): 


2-V2 

2 

0 

1 

0 

6  +  2V2 

7 

2-V2 

2 

0 

1 

73-34a/2 

-5  +  4a/2 

3-V2 

—  1  +  V2 

28 

4 

2 

2 

87-48V2 

-45  +  34a/2 

3-V2 

“1  +  ^2" 

28 

28 

2 

2 

(16) 


This  method  was  designed  to  have  a  stability  matrix  with  eigenvalues  zero  and  R(z),  where 


R(z)  = 


z{-j2-l) 

Z2(l-j2)-z(2--j2)  +  l’ 


the  same  stability  function  as  a  2-stage  Singly-Diagonally  Implicit  Runge-Kutta  (SDIRK) 
method  of  order  2  which  is  known  to  be  A-stable.  Therefore  this  method  (Equation  16)  is 
then  itself  A-stable. 

Type  3  (explicit,  parallel): 


'A  U~ 
B  V 


'  0  0 

0  0 
_3  _3 

8  8 
_7  9 

8  8 


1  0 

0  1 

2  1 
4  4 
3  7 
'4  4 


(17) 


This  method  has  a  stability  polynomial 


p(w,z)  =  w2  ~^l  +  ^w-j(3z  +  l), 

which  yields  a  stability  interval  along  the  real  axis  of  B4 
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Type  4  (implicit,  parallel): 


A  U 
B  V 


3-V3 

2 

0 

18-11^3 

4 

22-13^3 


0 

3  —  -s/3 


0 


0 


-12  +  7^3  3-2V3  -1  +  2-V3 


4 

-12  +  9-V3 


2 

3-2a/3 


2 

-1  +  2V3 


(IB) 


This  method  has  the  stability  polynomial 


p{w,z)  =  w7 


2  -  2V3  +  (8a/3  -  12)z  +  (15  -  9S)z7 
-2  +  (6  -  2V3)z  +  (3j3-6)?W  +  Z  (-2  +(6  -  2-75);  +  (3-73-  6)z2f 


2  +  (73-3)z 


and  is  found,  using  the  Schur  criterion,  to  be  A-stable. 

We  may  note  by  examining  Equation  2  that  for  the  serial  DIMSIM  types,  the  first 
internal  stage  must  be  calculated  before  the  second  stage,  etc.  For  the  parallel  types,  on  the 
other  hand,  the  calculations  of  the  internal  stages  are  independent  of  each  other  and  could 
conceivably  be  performed  on  separate  processors  of  a  parallel  computer.  This  is  an 
example  of  “parallelism  across  the  method,”  described  by  Gear  in  an  early  report 
(Reference  13)  in  which  this  form  of  parallelism  was  distinguished  from  “parallelism 
across  the  system”  in  which  different  equations  or  subsystems  might  be  handled  with 
different  processors.  Waveform  relaxation  (see,  for  example,  Reference  14  and  Reference 
15)  is  an  example  of  the  use  of  this  latter  form  of  parallelism  and  it  is  possible  that  both 
types  of  parallelism  could  be  combined  in  solving  a  single  problem.  It  is  evident  that  the 
parallelism  across  the  method  indicated  here  for  DIMSIMs  would  only  enable  effective  use 
of  a  number  of  processors  equal  to  the  number  of  internal  stages  while  parallelism  across 
the  system  could  employ  many  processors  in  solving  a  large  system.  It  may  also  be  noted 
that  die  parallel  types  have  fewer  parameters  and  it  would  then  be  expected  to  be  more 
difficult  to  find  methods  with  desirable  stability  properties  and  other  characteristics. 

The  type  3  methods  are  interesting  in  that  they  do  not  call  for  the  calculation  of  internal 
stages,  the  external  stages  are  the  same  as  the  internal  stages.  Also  since  A  =  0,  the  first 
equation  of  Equation  1 1  implies  that  the  external  stages  are  the  first  p  +  1  terms  of  the 

Taylor  expansions  of  y  at  the  stage  points  tn.,  +  C;h,  since  ecz  =  w{z)  +  o{zp+xy  The 

methods  then  use  a  Taylor  series  starting  method  and  time  marching  is  carried  out  using  the 
equation 
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(19) 


For  example,  with  p  =  2,  c,  =  0,  c2  =  1,  we  find  (using  methods  described  in  the  next 
section)  that  we  may  also  derive  the  method 


'A  U~ 
B  V 


0 

0 

0 

-1 

2 


0  1  O' 

0  0  1 

0  0  1 

1  0  1 

2 


(20) 


Upon  examining  the  time  marching  process  of  Equation  19  for  this  particular  method,  we 
see  that  it  is  equivalent  to  the  the  familiar  second  order  Adams-Bashforth  method, 


yn=yn-\+Wfn-\-fn-2)- 


For  Type  3  methods  the  stage  values  approximate  the  solution  at  the  stage  points,  since 
external  stages  and  internal  stages  are  equal  here  and  the  stage  order  condition  requires  that 
the  solution  be  approximated  by  the  internal  stages  at  the  stage  points.  These  are  each 
calculated  using  an  explicit  linear  multistep  method  depending  on  previous  values  only 
from  the  last  interval,  with  the  formula  for  the  ith  component  skipping  the  i  -  1  last 
calculated  solution  points.  A  D1MS1M  uses  an  interpolant  then  based  on  solution  and 
derivative  values  from  the  stage  points  to  produce  final  approximations  at  arbitrary  points. 
Nevertheless  the  Type  3  time  marching  process  can  be  seen  as  an  example  of  an  explicit 
cyclic  linear  multistep  method  (Reference  16)  when  stage  points  are  evenly  spaced,  and 
only  a  slight  modification  would  seem  necessary  for  irregular  spacing.  For  example,  the 
ith  element  of  the  external  stage  vector  approximates  the  solution  at  tn_,  +  c;h  and  is  given  by 


cjh,yf  1])+  l\ 


and  for  constant  step  size  and  uniform  spacing  of  c  from  0  to  1  (c-  =  (j-l)/(p-l))  and 
recognizing  that  we  really  have  a  solution  at  each  of  the  stage  points,  this  becomes  for  the 
ith  stage  of  the  nth  step,  using  the  simplified  notation  yk  to  the  refer  to  the  kth  stage, 
numbering  from  the  first  stage  of  the  first  step, 


=  h  £  yfo + (/>(»  -  2)  ■ +  y>./("‘2W )  +  X  vty 

7=1  7=1 


p{n-2)+j 
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and  this  is  in  the  familiar  explicit  linear  multistep  class  of  formulas.  Some  interesting 
features  of  this  family  when  seen  in  this  form  include  the  use  of  the  same  combination  of 
solution  values  for  each  element  within  a  cycle  (vTy[n‘n)  and  the  omission  of  several 
previous  values  and  derivatives  in  the  calculation  at  most  of  the  stage  points.  Clearly  this 
solution  process  can  utilize  separate  parallel  CPUs  for  the  p  derivative  evaluations  required 
at  each  step. 


APPROXIMATING  THE  SOLUTION  AND  THE  NORDSIECK  VECTOR 

Nordsieck  techniques  are  used  extensively  in  the  development  and  implementation  of 
DIMSIMs.  These  use  a  vector  of  derivatives  scaled  with  the  step  size  and  usually  also  with 
a  factor  of  l/(j  - 1)!  where  j  is  the  component  number  of  the  vector.  Here  we  omit  the  exra 
factor,  which  can  be  readily  provided  through  multiplication  by  a  constant  diagonal  matrix 
Q  =  diag(l/0!,  1/1!,  1/2!,  ...,  1/p!),  and  use  the  term  “Nordsieck  vector”  to  refer  to  a 
closely  related  vector  that  frequently  appears  here.  We  define  the  Nordsieck  vector  of 
length  p  +  1  as 


><*«) 

hy  {xn) 

hpy(p\xn) 


(21) 


Nordsieck  (References  17  and  18)  described  a  family  of  linear  multistep  methods  using  the 
vector  Qy(xn)  which  provided  an  especially  convenient  approach  to  step-size  changes. 

Nordsieck  noted  that  solutions  of  ordinary  differential  equations  could  be  reduced  to 
finding  interpolation  polynomials  to  use  in  representing  the  solution,  and  his  vector 
provided  a  readily  scalable  representation.  These  techniques  have  frequently  been 
incorporated  in  implementing  linear  multistep  methods  and  Butcher  and  Jackiewicz 
(Reference  3)  have  shown  how  they  may  be  utilized  effectively  with  DIMSIMs  as  well. 
Trocogna  (Reference  19)  has  since  then  used  this  approach  in  implementing  two-step 
Runge-Kutta  methods  (Reference  20). 

Two  matrices  are  defined  which  are  used  to  relate  the  Nordsieck  vector  to  the  internal 
and  external  stages  of  the  method.  Let  be  a  vector  with  kth  component 

Then  we  can  find  matrices  B  and  V  such  that 


y(x„)  =  +  Vy[n~l]  +  0(hp+l ). 


(22) 
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These  matrices  can  be  calculated  using  the  following  theorem,  announced  in  Butcher 
(Reference  1)  and  proven  rigorously  in  Butcher  and  Jackiewicz  (Reference  3).  We  first 
define  z  as  a  vector  of  length  p  +  1, 


1 


z  = 


z 


p 


(23) 


Theorem  3  (Butcher  and  Jackiewicz,  Reference  3):  Assume  that  the  method  in  Equation 
2  has  order  p  and  stage  order  q  =  p  or  q  =  p  -  1 .  Then  the  approximations  in  Equation  22 

are  correct  to  o{hp+l^  if  and  only  if 


ezz  -  zBecz  +  Vw(z )  + 


(24) 


Proof:  The  proof  was  developed  by  Butcher  and  Jackiewicz  and  is  reproduced  here  for  the 
sake  of  completeness.  Define  a  matrix  T  such  that 


1 


-1 


1 

2 


0  1 


-1 


(-l)p 

p! 


nr1 

(P7 1)!  ' 


0  0  0  •••  1 


Taylor  expansions,  the  stage  order  and  the  problem  definition  in  Equation  1  are  used  to 
obtain  the  relationships 


AXn) 

y(*n-l) 

hy'(xn) 

hpy{p\xn)_ 

0[hp+x), 


and 


17 


NAWCWPNS  TP  8340 


■«w 

hF(Y2) 

cq 

X^-i) 
hy'{xn- 1) 

hFiYp). 

L 

+  0(hq+1 ). 


We  also  use  Equation  3  to  write 


ii 

R 

N* 

if 

yfo-l) 

1 

'TT 

a 

i 

i 

+  0[hp+l). 


We  now  substitute  these  equations  into  Equation  22  to  obtain 

Be1'1  a 


P 

I 

z= 0 


hy%n-iW  =  +  f  W’Ov 

i=l  V*  i=0 


Since  q  =  porq  =  p-  lwe  can  combine  terms  to  obtain 

(Va0  -  tQ)y(xn_l)+  I  Vaf  +  — —  -  tt  y[l\xn_x)ti 

v  '  /=iv  / 

Equating  terms  of  the  same  power  in  h  we  obtain 
t0  =  Va0, 
and 

~  ck~l  ~ 

tk  =  B ^  _  i)!  +  k  — 


+o[hp+l)  +  o[hq+1). 


0{hp+l). 
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Multiplying  the  kth  equation  by  zk  and  adding  yields 
p  _  p  J-l-i-l  _  p 

ltizl=zBZ—^—  +  Vlaizl, 

i= 0  i= 1  (.*  -  1)!  i= 0 

which  may  be  seen  to  equivalent  to  Equation  29  upon  expansion  of  the  exponentials  in 
Taylor  series.  ■ 

Noting  that 

y[n]  =  Wy(xn)  +  0(hp+1) 

=  WBhF(Yln])  +  WVy[n~n  +  0(hp+l) 

=  BhF(Y[n])  +  Vy[n~l\ 


and  comparing  the  corresponding  terms  we  obtain  the  following  so-called  compatibility 
conditions 


WB  =  B, 
WV  =  V. 


(25) 


IMPLEMENTING  DIMSIMS 


IMPLEMENTATION  ISSUES 

In  order  for  an  ODE  solution  method  to  be  useful,  certain  capabilities  must  be  provided. 
A  numerical  method  will  most  effectively  be  applied  using  adaptive  step-size  selection 
based  on  an  error  tolerance.  Thus  neither  too  much  work  is  done  due  to  steps  that  are  too 
short,  nor  is  required  accuracy  sacrificed  by  using  steps  that  are  too  long.  This  requires  an 
ability  to  change  step  size,  and  also  the  ability  to  estimate  error  and  suitable  step  length. 
Order  changing  is  also  desirable  for  maximum  efficiency,  but  techniques  for  accomplishing 
this  with  DIMSIMs  have  not  yet  been  developed.  Furthermore,  although  typically  the  ODE 
solution  is  desired  at  certain  output  points,  often  evenly  spaced,  the  ideal  combination  of 
efficiency  and  accuracy  calls  for  integration  steps  to  be  as  long  as  accuracy  will  allow. 
Thus  interpolation  should  be  used  to  obtain  output  values.  Some  sort  of  starting  technique 
is  required  for  the  first  step.  Finally,  for  methods  for  which  c,  =  0  and  cp  =  1,  a  significant 
savings  in  work  is  possible  by  using  a  new  approach  that  will  be  described  for  evaluation 
of  the  first  internal  stage. 
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CHANGING  STEP  SIZE  WITH  DIMSIMS 

We  must  introduce  some  notation  to  discuss  the  rescaling  that  becomes  necessary  to 
continue  the  calculation  with  a  modified  step  size.  We  begin  by  calculating  y[1]  using  Ym 
and  h,  =  t,-t0,  and  we  assume  now  that  y[n'rt  has  been  calculated  usingY[n‘irand  hn.j  =  tn_,  - 
t„-2-  We  now  calculate  a  y[n]  corresponding  to  tn  using  Y[n]  and  step  size  hn  =  tn  -  tn_j.  We 
denote 


fo-i)  = 


(26) 


We  also  define  the  diagonal  matrix 


D  =  diag(l,S,S2,...,SP),  (28) 

where  8  =  dn=— — .  Furthermore  we  distinguish  y^~^and  y^”-1^  where  y^"-1]  is 
hn-l 

rescaled  to  reflect  a  next  step  of  hn.  This  rescaling  process  will  be  described  below.  We 
desire  to  have  approximate  in  the  same  way  that  yt"-1^  approximates 

Wy(tn_ 1).  This  leads  to  the  relations 
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j*"”1'  =  Wy(<«-\)+o{K-t)- 

= K-iBF[y{n~'])+ V51"-21 + ohd),  v[«l = yl°l, 
=  WDy(t„_ i)  =  J1"-11  +  0(hP"). 


This  indicates,  after  using  the  expression  for  the  Nordsieck  vector  at  tn_,,  that  we  can 
rescale  by  using  the  formula 


y[n~1]  =  hn_xWDBF{Y{n~X])  +  WDVy[n~2]. 


(29) 


We  now  have  a  modified  numerical  process,  as  follows: 


(30) 


yM  =  /i„wMf(fW)  +  WDVy^n~l  l 


We  note  (Reference  1)  that  zero-stability  for  a  step-changing  solver  will  then  be  determined 

by  the  eigenvalues  of  the  matrix  WDV .  Thus  any  free  parameters  in  determining  B  and  V 
might  well  be  used  to  ensure  good  stability  for  step-size  changing.  That  this  is  always 
possible  will  now  be  shown. 

Computation  of  the  rescaling  matrices  is  simplified,  first  of  all,  with  all  type  1 
DIMSIMs  with  c,  =  0.  The  equations  to  be  satisfied  include  Equation  25  (the  compatibility 
condition),  Equation  24  and  the  condition  that  step-size  changing  zero  stability  be 

nonrestrictive,  met  where  the  matrix  WD  V  has  only  the  one  nonzero  eigenvector  e  with  the 
associated  eigenvalue  of  1.  We  are  able  to  meet  both  of  these  conditions  for  a  choice  of  V 

such  that  the  first  row  is  vT  and  all  the  other  rows  are  0,  and  for  a  choice  of  B  such  that  the 
first  row  is  the  same  as  the  first  row  of  B  and  the  other  elements  are  computed  uniquely 

from  a  linear  system.  We  quickly  note  that  for  V  of  this  form,  DV=  V,  and  WV=V  by 
the  compatibility  condition,  and  V  is  defined  to  be  equal  to  evT.  This  indeed  has  the  one 
nonzero  eigenvector  e  with  associated  eigenvalue  of  1.  We  now  examine  the  condition  24 
for  the  first  rows  of  the  two  rescaling  matrices. 
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ez  =  zBxecz  +  V]w(z)  +  o[zp+l ). 


We  compare  this  with  the  similar  condition  in  Equation  1 1  for  the  first  rows  of  B  and  V, 
and  noting  on  the  left  side  the  condition  c,  =  0,  which  makes  Wj(z)  =  1,  we  obtain 


=  ZB1ecz  +  VMz)  +  o(zp+1). 


Clearly  the  first  rows  satisfy  the  same  conditions,  and  hence  the  first  rows  of  B  and  V  can 
be  taken  the  same  as  the  first  rows  of  B  and  V,  respectively.  We  now  consider  the 

conditions  on  the  remaining  rows.  If  we  take  the  remaining  rows  of  V  to  be  0,  we  then 

have  the  conditions  for  rows  of  B  as  follows: 


ezzk~X  =zBkecz+0{zp+l\  k  =  2,...,p  +  l. 


We  have  p  unknowns  in  each  equation.  But  if  we  examine  the  number  of  conditions 
required  to  equate  polynomial  terms  of  degree  0  through  p,  we  find  that  there  will  be  no 
degree  0  condition,  but  degrees  1  through  p  will  always  appear.  Thus  we  have  p  linear 

equations  in  p  unknowns  for  each  remaining  row  of  B.  Also,  the  matrices  for  each 
equation  are  the  same  but  the  right  hand  sides  will  always  be  different.  And  so  long  as  the 
method  is  nonconfluent  (that  is,  c,-  ^  cj  for  i  ^  j),  there  should  be  no  problem  in  solving 

these  systems,  as  may  be  observed  from  the  formation  of  the  matrix  elements  from  the 
Taylor  series.  This  is  easily  extended  for  both  implicit  and  explicit  DIMSIMs  when  c1  is 

not  0.  Taking  V  as  before,  we  find  the  situation  with  B  for  rows  after  the  first  to  be 
unchanged.  The  relationship  of  the  equation  for  the  first  row  to  the  equation  for  the  first 
row  of  B  is  modified,  however,  since  Wj(z)  is  no  longer  1.  Recognizing  that  the  the  zero 
degree  term  in  each  component  of  w(z)  is  1,  we  find  that  the  equation 


ez  -  vTw(z )  =  zB\ecz  +  o{zp+X ) 


has  a  zero  term  of  degree  0,  since  the  leading  1  in  the  Taylor  expansion  of  e2  is  cancelled  by 
the  leading  1  from  the  term  vTe.  Thus  we  again  have  p  equations  in  p  unknowns.  Note 

that  the  compatibility  condition  for  B  must  automatically  be  satisfied  for  this  form  of  V 
since  multiplying  Equation  24  by  W  and  setting  6  =  1  produces  the  same  equation  that  B 
must  satisfy,  the  second  relation  of  Equation  11.  We  summarize  these  observations  in  the 

following  theorem  which  will  be  helpful  in  calculating  the  rescaling  matrices  B  and  V. 
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Theorem  4:  For  nonconfluent  DIMSIMs,  we  have 

i)  the  rescaling  matrix  V  may  be  chosen  to  consist  of  a  first  row  vT  and  all  other  rows 
equal  to  0; 

ii)  if  Cj=0,  the  rescaling  matrix  B  will  then  have  a  first  row  that  is  identical  to  the  first 
row  of  B, 


ezzk~l  =  zBkecz  +  o(zp+l),  k  =  2,...,p  +  U 


iii)  if  Ci  &  0,  the  first  row  of  B  must  satisfy  the  linear  system 
ez  - vTw(z )  =  zB\ecz  +  o{zp+^\ 


iv)  and  for  all  choices  of  cp  subsequent  rows  must  each  satisfy  the  linear  system 
resulting  from  Taylor  expansion  of 


ezzk~l  =  zBkecz  +  o(zp+1),  k  =  2,...,p  +  l. 

We  see  that  the  coefficients  of  B  for  rows  after  the  first  are  thus  dependent  only  on  the 
choice  of  stage  points.  It  is  evident  that  the  higher  derivatives  of  the  Nordsieck  vector  are 
obtained  in  this  approach  by  taking  advantage  of  the  high  stage  order.  Applying  the 
derivative  function  to  the  internal  stages  yields  high  order  approximations  to  die  solution 
derivatives  at  the  stage  points  which  yield  linear  systems  for  the  higher  derivatives  through 
Taylor  expansion . 

Zero  stability  is  obviously  a  rather  mimimal  condition.  More  careful  step-size  change 
stability  requires  examination  of  the  stability  matrix  for  this  process.  This  may  be  written 
as 


M(S)  =  zWD(S)B{I  -  zA)~l  +  WD(8)V  =  zWD(S)B(I  -  zA)-1  +  V  (31) 


for  the  customary  form  of  V  given  by  theorem  4.  But  it  is  not  true  that  eigenvalues  of  this 
matrix  will  determine  the  growth  of  the  external  stage  vector,  since  it  is  nonsymmetric  and 
typically  varies  with  each  step,  and  this  becomes  a  very  difficult  problem.  However, 
lacking  some  better  criterion,  in  a  heuristic  approach  the  eigenvalues  of  this  matrix  may  be 
examined  by  calculating  a  sort  of  “pseudo-stability  region”  to  provide  some  indication  of 
the  effect  on  stability  of  step-size  changing  with  a  view  to  aiding  the  determination  for  a 
solver  of  the  bounds  to  set  on  step-size  changes.  This  approach,  pursued  by  Enenkel 
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(Reference  21)  in  his  study  of  related  general  linear  methods,  yielded  very  restrictive 
results. 


BUTCHER- JACKIEWICZ-TYPE  INTERPOLATION 

Butcher  and  Jackiewicz  (Reference  3)  proposed  continuous  interpolants  of  uniform 
order  p  of  the  form 


ri(tn- 1 +QK)  =  KM^)F(Y[n]) + ro(9)y[n~1]- 


(32) 


Here  we  define  p0(9)  =  [pQl(6),f302(d),...,f30s(d)]  and  yo(0)  =  [7oi(0)>7o2(0)>->7or(0)]’ 
where  the  components  are  polynomials  of  degree  p  (or  lower  if  certain  coefficients  become 
set  to  zero).  For  a  Nordsieck  interpolant  we  may  write  the  interpolant  in  this  form  as  well, 

and  with  the  customary  form  for  V,  Yq{9)  =  v  .  We  further  note  that  in  order  for 
compatibility  with  the  equation  for  the  first  component  of  the  Nordsieck  vector,  /30(1)  and 

70(1)  must  be  equal  to  the  first  rows  of  B  and  V,  respectively.  This  interpolant 
compatibility  condition  will  be  shown  to  be  incorporated  within  the  order  condition 
provided  in  the  following  theorem  derived  by  Butcher  and  Jackiewicz  (Reference  3)  in 
which  the  case  q  =  p  -  1,  though  not  of  particular  interest  here,  is  also  included. 

Theorem  5  (Butcher  and  Jackiewicz):  If  a  DIMSIM  has  order  p  and  stage  order  q  =  p  or 
q  =  p  -  1 ,  then  rj  approximates  y  with  uniform  order  p  if  and  only  if 


zpQ(d)ecz  +  y0(6)w(z)  =  +  0{zp+l),  6  €  (0,1], 


(33) 


and  w(z)  is  as  defined  above.  Moreover,  the  interpolant  r\  is  continuous  on  the  whole 
interval  of  integration  if  and  only  if 


4(0)  =  0, 

Yo(0)WDB  =  /?0(1)>  (34) 

7o(0)WZ>F  =  7o(1). 


Proof:  Expanding  y(tn),  and  y’(tn+]+ch)  =hF(Y[nl)+0(hp+I)  in  Taylor  series  about  tn_j 
and  assuming  that  the  interpolant  approximates  y  to  uniform  order  p,  we  may  write 
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f  = AM  X^f-AV^-i) + r»(«)  i + o(*r') + o(*r2). 

*=o  Jfe!  Jt=i  (A  —  1)!  *= o 


jfc-1 


where  we  denote  the  (i-l)th  column  of  W  as  a,.  Since,  as  previously  noted,  the  first 
column  of  W,  a0,  is  e,  and  q  =  porq  =  p-  l,we  can  write 


(yo(0)e-i)y(xn-i)+  I 

k= l 


f 


M)£ 


k\ 


^•+ro(e)«t  ~\kym{*n-i)= o(h^). 


We  then  generate  p  +  1  equations  by  setting  coefficients  of  powers  0  to  p  of  hn  to  0: 


7o(0)e  =  l, 

j3o(0)(F^+7o(eK  =  ¥’  l-k~p- 


We  may  now  multiply  each  equation  by  an  appropriate  power  of  z  and  sum,  obtaining  an 
equation  equivalent  to  the  result  we  seek: 


P  ck-xzk~l 


k= 1  (*-!)! 


+  y0{6)Zakzk 


k=0 


p  6^ 
k= 0  fc! 


These  steps  are  reversible. 

For  continuity  at  node  point  xn  we  must  have  7](xn  -)  =  ri(xn  +).  The  equation 
representing  this  relationship  reads: 


MoW/fi’1"1) + roO)}1”'11 = A+iAfo)/^”*11) + ro(o)yw 

= + y^iy/DBf{^) + wdv?1”'11). 


Comparing  coefficients  of  f(Y[n+1])  we  see  that  P0(0)=0.  Similarly,  comparison  of 
coefficients  of  f(Y[nl)and  y[n'1]  yields  y0(0)WDB  =  /30(1)  and  y0(0)WDV  =  7oC0> 
respectively.  ■ 

The  following  corollary  provides  simplification  of  the  conditions  of  this  theorem. 
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Corollary:  If  a  DIMSIM  has  order  p  and  stage  order  q=porq  =  p-  1,  then  77 
approximates  y  with  uniform  order  p  if  and  only  if 


zPo(B)ea  +ro(SMz)  =  e*  +o(z’,+l),  9e  (0,1], 


(35) 


and  w(z)  is  as  defined  above.  Moreover,  the  interpolant  77  is  continuous  on  the  whole 
interval  of  integration  if  and  only  if 


A(0)  =  0. 

7o(O)£  =  /30(l),  (36) 

7o(0)V  =  7o(l)' 


Furthermore,  if  these  conditions  are  met,  the  interpolant  compatibility  conditions  that  /?0(1) 

and  y0(l)  must  be  equal  to  the  first  rows  of  B  and  V,  respectively,  will  automatically  be 
satisfied. 

Proof:  Looking  again  at  the  order  condition  in  Equation  33,  we  find  that  because  of 
continuity,  we  can  now  consider  the  case  where  6  is  set  to  0.  Then  since  /J0(0)  =  0,  and 

e0z  =  l,  we  have  70(0)w(z)  =  1.  Now  w(z)  =  Wz.  We  can  consider  the  various 
polynomial  terms  separately.  Let  ex  e2,  etc.,  be  the  unit  vectors  with  Is  in  the  appropriate 
position.  Then  y0(O)W£j  =  l,  and  YQ(<S)Wek  =  0  for  k  >  1.  Thus 

70(0 )WD  =  [l  0  •  •  •  0]  regardless  of  the  value  of  8.  Thus,  from  the  order  condition  we 

have:  Yq(0)WDB-/5q  (the  first  row  of  B)  and  7q(0)1VDV  =  y0  (the  first  row  of  V). 
Note  that  the  interpolation  coefficients  will  never  depend  on  8.  Also  note  this  eliminates 
the  interpolant  compatibility  conditions  as  separate  criteria;  since  there  is  no  dependence  on 

its  value  we  may  arbitrarily  set  8  =  1,  in  which  case  D  =  I.  Then  WDB  =  WIB  =  WB  =  B, 
and  WDV  =  WIV  =  WV=V.  ■ 

Deriving  an  interpolant  is  then  a  matter  of  finding  coefficients  to  satisfy  these 
relationships.  In  this  report  all  interpolants  of  the  form  in  Equation  32  and  satisfying  the 
conditions  of  Theorem  5  will  be  termed  Butcher-Jackiewicz-type  interpolants,  while  the 
interpolants  proposed  below  by  the  author  of  this  report  will  be  termed  Nordsieck  or 
continuous  Nordsieck  interpolants. 

Although  examples  of  continuous  interpolants  of  maximal  order  will  be  derived  and 
utilized  in  the  Implementing  DIMSIMS  section,  we  must  note  here  that  a  continuous 
interpolant  of  maximal  order  in  this  form  does  not  exist  for  all  DIMSIMs.  For  example,  we 
consider  the  type  2  example  of  Equation  16,  shown  again  below: 
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2-V2 

2 

0 

1 

0 

6  +  2^2 

7 

2-^2 

2 

0 

1 

73-34^2 

-5  +  4a/2 

3-V2 

-1  +  V2 

28 

4 

2 

2 

87-48^2 

-45  +  34^2 

3-^2 

-1  +  V2 

28 

28 

2 

2 

with  c=[0,l]T.  For  this  method  the  matrices  B  and  V  were  calculated  by  Butcher 
(Reference  1)  to  be 


'-26  +  41V2 

62-37^2" 

28 

-48  +  51V2 

OO 

64-33^2 

28 

-20+15V2 

90 

20-9^2 

14 

14 

'-12+1W2 

26-1  lV2~ 

14 

-30  +  3a/2 

14 

30-3^2 

14 

12  +  3^2 

14 

12  +  3^2 

7 

7 

(Note:  A  separate  calculation  using  the  customary  form  for  V  yielded  the  same  result.) 

We  now  apply  the  compatibility  and  continuity  conditions  (with  the  help  of  Mathematica) 
and  show  that  an  interpolant  of  the  form  described  here  cannot  satisfy  these.  We  look  for 
vectors 


A)(^)  “  [AhO  +  + A)12^2’A)20  +  A)210  + A)22^2]’ 

7o(^)  =  [7oio  +7oh^  +  7oi2^2’7o20  +7o2i^  +  7o22^2]- 
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We  immediately  note  from  the  first  compatibility  condition  that  P010  and  (3020  must  be  0. 
We  use  the  compatibility  conditions  to  eliminate  (3012,  p022,  y011,  and  y021.  We  then  use 
the  second  continuity  condition  to  eliminate  y012  and  y022.  But  then  the  third  continuity 
condition,  y0(O)V=yo(l),  becomes  equivalent  to 


1 142683156482  - 1429914988166^2 
87418383556 


-6579484927094450  -  4651660909947086^2 
87418383556 


=[0  0], 


which  contradicts  the  assumption  that  an  interpolant  of  the  given  form  exists. 

A  consequence  of  the  fact  that  interpolants  of  the  given  form  do  not  always  exist  is  that 
a  search  for  a  suitable  DIMSIM  scheme  should  incorporate  conditions  for  the  existence  of 
an  interpolant  of  this  desirable  form.  The  existence  of  a  suitable  interpolant  is  crucial  for 
the  implementation  of  DIMSIMs  in  a  waveform  relaxation  strategy,  and  for  those  methods 
for  which  this  form  does  not  exist,  the  Nordsieck  form  may  always  be  used.  A  subsequent 
report  will  compare  the  performances  of  alternative  interpolants  for  DIMSIM 
implementations  for  waveform  relaxation. 

We  may  readily  obtain  the  following  corollary  to  Theorem  5: 

Corollary:  If  a  continuous  interpolant  of  the  form  of  Equation  34  exists,  then  there  exists 
a  constant  vector  70(0)  such  that 

i)  y0(0)y^  =  yn  for  each  step  of  the  integration  process,  where  yn  is  the 
approximation  to  the  solution  at  tn  and  is  also  given  by  the  first  component  of  the  Nordsieck 
vector  calculated  using  Equation  22. 

ii)  70(0)W  =  [1  0  0], 

iii)  70(0)5  =  B{,  y0(0)V  =  V\  =  V\  =  vT,  where  the  subscript  1  indicates  the  first 
row  of  the  matrix  and  V  =  evT. 

Proof:  i)  Setting  0  to  0  in  Equation  34  and  using  continuity  at  tn_,  and  /J0(0)  =  0,  we 
obtain 


rt*n- 1)  =  7o(%["  1]  =  V(tn- 2  +  Vi)  =  7i(?«-i)  • 


We  then  note  that  we  may  simply  adjust  the  subscript, 
ii)  See  proof  of  Theorem  5. 
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iii)  Since  the  first  column  of  W  is  e,  from  ii),  y0(0)e  =  1.  Then 
y0(0)y  =  y0( 0)evT  =  vT .  Now  see  proof  of  Theorem  5.  ■ 


NORDSIECK  INTERPOLATION 

The  availability  of  the  Nordsieck  vector  provides  a  ready  interpolant.  For  any  DIMSIM 

we  may  calculate  coefficient  matrices  B  and  V  using  theorem  4,  which  will  then  yield 
Nordsieck  vectors  at  each  grid  point.  Interpolation  can  be  carried  out  backward  and 
extrapolation  forward  from  a  grid  point  using  the  Taylor  expansion  polynomial  of  degree  p. 
Since  the  Nordsieck  vector  components  are  all  locally  accurate  to  0(hp+1),  the  global 
accuracy  of  the  interpolant  is  then  0(hp),  the  same  as  the  accuracy  of  the  method  at  the  grid 
points.  This  is  provided  for  all  DIMSIMs,  including  the  Type  2  method  above  for  which 
the  Butcher- Jackiewicz-type  interpolant  does  not  exist.  (Note  that  the  availability  of 
extrapolation  accurate  extrapolation  forward  is  extremely  useful  in  that  it  provides  an 
excellent  internal  stage  predictor  for  implicit  methods.  This  will  be  developed  more 
extensively  in  a  report  in  preparation  on  implicit  DIMSIMs.)  If  we  have  calculated  an 
approximation  to  the  Nordsieck  vector  at  tD  as 


we  may  then  carry  out  Nordsieck  interpolatation  at  tn_ j  +  91^  using  the  Taylor  series 
formula  at  the  point  tn  of  the  form 


7](r„_i+0A„)  = 


1  -1 


p\ 


m-0)y{tn ). 


where 


'1 

0 

0  • 

•  0  ‘ 

0 

8 

0  • 

■  0 

0 

0 

82  ■ 

•  0 

0 

0 

0  • 

•  $P_ 

(37) 


and  where  1-0  is  used  because  interpolation  is  carried  out  to  the  left  of  the  grid  point,  at  tn- 
(l-0)hn.  Alternatively  interpolation  to  the  right  can  be  carried  out.  The  formula  changes 
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only  slightly,  with  0  used  instead  of  1-0  and  the  negative  signs  in  the  multiplying  vector 
are  eliminated. 

This  interpolant  has  two  drawbacks.  First,  since  a  fresh  DIMSIM  calculation  of  the 
Nordsieck  vector  is  carried  out  at  each  grid  point,  these  will  be  points  of  discontinuity. 
Secondly,  Taylor  series  decrease  in  accuracy  away  from  the  node  point.  Both  of  these 
problems  can  be  overcome,  or  at  least  mitigated,  by  increasing  the  Taylor  series  degree  to 
p  +  1.  The  extra  component  is  calculated  to  exactly  remove  the  discontinuity  at  the  far  grid 
point  and  will  yield  a  0  contribution  at  the  grid  point  where  the  Nordsieck  vector  is 
calculated.  For  example,  the  interpolant  above  would  be  modified  to  produce  a  continuous 
Nordsieck  interpolant  by  the  addition  of  the  term 


(Wn-i)-  1 


V 


(38) 


Since  this  term  should  be  quite  small  due  to  the  higher  power  of  hn,  the  overall  behavior 
of  the  interpolant  should  not  degrade,  and  in  fact  tests  described  below  in  the  chapter, 
“Developing  Explicit  DIMSIM  ODE  Solvers,”  indicate  that  it  performs  quite  well. 


DIMSIM  ERROR  ESTIMATES 

An  error  estimate  is  intended  to  approximate  the  local  error  in  the  solution  of  an  initial 
value  problem.  That  is,  if  a  DIMSIM  has  produced  a  solution  at  grid  point  tn.,of  yn.,,  then 
the  DIMSIM  solution  yn  at  tn  is  compared  to  the  exact  solution  y^t^.y^)  of  the  problem 


1/(0  = /MO) 

1  y^n-i) = yn-\  ’ 


and  the  local  error  at  tn  which  is  to  be  estimated  is  defined  to  be 


errn=yn-y(tn]tn_hyn_l). 


(40) 


Unlike  linear  multistep  and  Runge-Kutta  methods,  DIMSIMs  do  not  typically  produce 
an  approximation  to  the  solution  at  the  grid  points  in  the  time  marching  process.  The  first 
component  of  the  Nordsieck  vector  must  be  separately  computed  from  the  external  and 
internal  stages  when  an  approximation  to  the  solution  is  desired.  Furthermore,  the  external 
stages  are  neighboring  trajectories  that  are  averaged  in  a  special  way  to  produce  a  solution 
and  the  next  internal  and  external  stage  vectors.  Finally,  the  Nordsieck  vector  components 
deviate  essentially  independently  as  the  solution  process  proceeds.  Thus  after  the  first  step, 

the  kth  component  approximating  hky^k\tn)  does  not  in  fact  represent  to  0(hp+1)  the  scaled 
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kth  derivative  of  the  solution  of  the  ordinary  differential  equation  y’(t)  =  f(t,y)  passing 
through  the  point  y(tn)  given  by  the  first  component,  but  rather  to  0(hp).  Thus  error 
estimation  may  be  expected  to  proceed  differently  from  the  approach  in  the  older  families  of 
methods.  However,  the  work  of  estimating  local  error  is  still  to  determine  the  amount  by 
which  the  approximation  to  the  solution  at  the  end  of  a  step  deviates  from  the  local  solution 
of  the  initial  value  problem  y’(t)  =  f(t,y)  passing  through  the  approximation  to  the  solution 
obtained  at  the  end  of  the  previous  step.  And  this  has  proven  to  be  possible  with 
DIMSIMs. 

Although  error  estimates  for  other  cases  have  been  derived  (Reference  3),  we  will 
restrict  ourselves  to  the  case  p  =  q.  The  stage  order  condition  requires  that  if  the  solution  is 
sufficiently  smooth,  we  may  write 

%  =  +  c,h)  +  +  0[h"*2),  (41) 


where  the  middle  term  on  the  right  is  designated  the  principal  part  of  the  error  for  Yj.  We 
note  here  that  we  are  using  the  function  y  to  refer  to  the  exact  local  solution  described 
above.  Although  theoretical  investigation  of  the  effect  of  the  accumulation  of  global 
solution  error  on  the  validity  of  this  local  stage  order  condition  remains  to  be  carried  out, 
there  is  ample  experimental  evidence  to  support  the  conjecture  that  Equation  41  holds  true, 
and  that  global  errors  do  not  reduce  the  order  of  the  leading  stage  error  term.  We  assume 
that  this  is  the  case  in  the  following  development,  which  otherwise  closely  follows  the 
approach  developed  by  Butcher  and  Jackiewicz  (Reference  3),  in  which  possible  stage 
order  reduction  is  taken  into  account. 

We  note  that  the  following  general  definition  for  local  discretization  error  of  the  external 
stages  applies  to  all  DIMSIMs.  The  idea  is  to  identify  what  is  to  be  called  the  local 
discretization  error  of  the  external  stages  with  the  term  of  order  hp+1  in  the  difference 
between  the  exact  external  stages  at  tn  and  the  calculated  external  stages,  assuming  that  an 
exact  Nordsieck  vector  is  used  initially  at  tn.,  and  that  stage  and  order  conditions  are  met. 

Definition:  The  local  discretization  error  le;(tn)  of  the  ith  external  stage  yj^of  the 
method  1(2.2)  at  the  point  tn  is  given  by 


h(‘n)=  £  £vsaJt/)(r„_,)A*.(42) 

k= 0  k=\  '  '  j=lk=0 


where 


4nl=h£  aijf(tn.l+cjh,r}"])+  X  ,/>(«_, )A\  *  =  1.2.....S.  (43) 

j= 1  '  '  7=1  (=0 
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For  the  case  p  =  q  =  r  =  s,  U  =  I,  and  assuming  that  the  local  condition  for  stage  order 
holds,  we  make  these  substitutions  in  Equation  42  to  obtain  the  simplified  expression 


&=0  7=1  v  1 


.(44) 


In  any  case,  the  vector  of  values  le^tj  we  designate  as  the  local  discretization  error  of  the 
external  stages. 

Theorem  6  (Butcher  and  Jackiewicz,  Reference  3)  The  local  discretization  error  le(tn)  of 
the  external  stages  of  Equation  2  at  the  point  tn  is  given  by 


le(x„)  =  'Ppy(p*'\xn-l)hpH  +  0(*',+2),  (45) 


where 


(Pp  = 


P+'Up+l-k 

t\  k\ 


BcP_ 
pi  ‘ 


(46) 


Proof:  The  proof  is  given  for  more  general  choices  of  p,  q,  r,  and  s  in  Butcher  and 
Jackiewicz  (Reference  3)  and  is  reproduced  here  for  this  restricted  case  for  the  sake  of 
completeness.  We  may  use  Taylor  expansion  to  express 


/>(<„)  =  '  +  0(hr*2-k\k  =  0,1, 

1=0  p 


P+1. 


We  then  substitute  this  expression  into  Equation  44  and  obtain 

l  =  hi  +  cjh) + 

k= 0  1=0  "■  j= 1 

f  +‘e,(tn)  +  o(h’‘*2). 

j=lk=0 


32 


NAWCWPNS  TP  8340 


We  now  expand  y’(tn.,+Cjh)  about  tn.j  to  get 


p  p+\-k 


k= 0  1=0  11  7=1  k=l  \K  V* 


+it,^yy,-,y+ie,{t,)+o(h^). 

j= 1  k=0 


We  now  reorder  the  summation  on  the  left  side: 


p  p+l-k 


XX  K' 

k=0  1=0  l- 


P+ 1 


*=o/=0  IV -/cj!  /=i  /! 

= xfx^H>‘ +fPJfiyir'%-t)h’M- 

k=0\l=0  1 '  J  1= 1  L 


Then,  interchanging  summation  orders  on  the  right,  we  obtain 


(  k 


X'ai,Pj±±(P+ 1)/,  \>p+l 


x  x^  n'.-,K +x^<'*,,(v,> 

jt=oV/=o  y  /=1  Ll 

= ft  £*>(0** +XX  w”('.-,)+ +  ofA'*1). 

A=1  j= l  ”  if!  *=0  7=1 


We  now  combine  tenns  and  rearrange  to  find  that 


P  (  k  M n  p  P  'N  ljk 

x  x^-x^vr- -x«v *1. 

jt=ov/=o  L  7=i  7=i  y  /c- 


+ 


x^  -  x^4 

\^/=l  "  ;=1  P-  J 
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The  order  conditions  ensure  that  terms  of  order  less  than  p  +  1  in  h  must  vanish,  and  so  we 
have  the  result  that 


P£  ai,p+l-l 
1=1 


l\ 


P 

"I 

7=1 


bjjC?  A 
p]- 


M+l\tn_x)hp+l  +0[hp+l). 


The  connection  between  errors  in  the  external  stages  and  errors  in  the  solution  values  is 
not  immediately  obvious.  However,  the  value  of  wlelxj  takes  on  a  special  significance. 
Here  vT  is  both  the  row  vector  (all  identical)  of  V  and  the  left  eigenvector  associated  with 
eigenvalue  1  of  V.  This  will  be  called  the  principal  part  of  the  local  discretization  error. 
Albrecht  (Reference  22)  showed  that  for  a  broad  class  of  methods,  this  is  the  quantity  that 
should  be  controlled  to  maintain  accurate  integration,  and  as  will  be  demonstrated  in  the 
Techniques  for  Obtaining  Values  for  DIMSIMs  section,  this  choice  leads  to  very 
satisfactory  results. 

The  error  estimate  for  fixed  step  sizes  may  be  found,  as  demonstrated  by  Butcher  and 
Jackiewicz  (Reference  3),  in  the  form 


vT(ppy{p+1\tn_l)hp+1=hf3TF(Y[n])  +  rTy[n~1]+o(hp+2), 


(47) 


and  for  variable  step  sizes  only  the  minor  modifications  shown  below  are  needed.  The 
error  estimate  is  then  a  linear  combination  of  terms  aready  computed  and  thus  is  essentially 
free  of  computational  cost.  The  (3  and  y  vectors  are  determined  by  the  method  and  may  be 
different  for  initial  steps,  constant  step  sizes,  and  varying  step  sizes  (involving  8  but 
reducing  to  the  constant  step  size  formula  for  8=1).  They  may  be  computed  as  follows. 

Define  a  matrix 


G  = 


(48) 


and  a  matrix 
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1  -1  i 


0  1  -1 


0  0  0 


(-1) 


P+1 


(p  + 1)! 

HI 

p} 

l 


(49) 


Then  the  following  theorem  applies: 


Theorem  7  (Butcher  and  Jackiewicz,  Reference  3):  If  the  fixed  step  method  in  Equation  2 
has  stage  order  q  equal  to  the  order  p  and  V  is  a  rank  one  matrix  such  that  Ve  =  e,  then 
vTle(xn)  can  be  estimated  by  the  formula 


vTle(xn)  =  hpTF[Y[n])  +  yTy[n~1]  +  o(hp+2). 


(50) 


where,  for  every  step  after  the  first,  the  vectors  P  and  y  satisfy  the  system  of  equations 


yTe  =  0 

+  yTBGti=  0,  i  =  l,2, ...,p  . 

(i-l)! 

pT£L  +  rTBGl  T 

pl 


(51) 


Proof:  The  following  is  based  on  an  earlier  proof  of  Butcher  and  Jackiewicz  (Reference  3) 
and  is  included  here  for  completeness  and  to  indicate  where  their  proof  requires  that  the 
step  number  n  be  greater  than  1 . 

The  first  equation  of  Equation  51  is  the  equation  resulting  from  terms  of  0(1)  in 
Equation  50,  since,  as  we  have  already  noted,  the  first  column  of  W  is  e,  and  so  each 
component  of  ytn'1]  is  y(xn4)  +  O(h).  Since  each  column  of  a  rank  1  matrix  is  proportional 
to  e  and  V  is  chosen  to  be  rank  1,  we  must  then  also  have  yTV  =  0.  We  now  use  Taylor 
expansion,  the  stage  order,  and  the  initial  value  problem  Equation  1  to  write 
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hy'{xn- 1) 


+  o[hp+ 2) 


hP"y^'\x„A 


We  may  also  write 


+  o(hp+1). 


hP+'yir+'XxA 


but  we  note  that  this  only  makes  sense  for  n  >  1,  that  is,  for  steps  after  the  first.  We  may 
also  write  the  Taylor  series  relationship 


1 

<N 

1 

y(xn- l) 

hy'(xn_2) 

=  f 

•• 

0(hp+2), 


which  of  course  only  applies  to  the  first  step  if  the  domain  of  validity  of  the  differential 
Equation  1  extends  sufficiently  far  to  the  left  of  x0.  Finally  we  use  the  second  equation  of 
the  method  Equation  2  to  write 


+  Vy[n~2\ 


which  clearly  is  meaningless  for  the  first  step.  We  substitute  these  expressions  into 
Equation  50,  and  using  yTV  =  0  and  Theorem  6  we  obtain 


vr<^+1y(;?+1)(xn_i)  =  (pTG  +  yTBGf) 


y{*n- 1) 

l) 

hp+ly{p+%n.  i). 


+  0{hp+1).  (52) 
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But  we  have 


pTG  +  yTBGf  =  0  ,pTe  +  yTBGih...,^-^-  +  yTBGt +1 

PI 


Substituting  into  Equation  52  and  equating  terms  of  the  same  degree  in  h,  we  obtain  the 
conditions  of  Equation  51.  ■ 

For  the  first  step,  the  accuracy  of  the  approximation  used  in  generating  the  starting 
vector  y[0]  becomes  important.  Although  the  minimum  requirement  is  that 


^]=laiky{t)(h)hk  +  o{h-*'), 
i=0 


it  is  possible  to  obtain  more  accurate  estimates  for  the  higher  derivatives  of  the  solution  so 
that  the  accuracy  becomes  0(hp+2).  The  fixed-step  error  formulas  derived  by  Butcher  and 
Jackiewicz  may  then  be  replaced  by  formulas  of  similar  form  but  with  different 
coefficients. 

Theorem  8  (Initial  Step  Error  Estimate):  If  the  solution  y(x)  to  the  problem  in 
Equation  1  is  sufficiently  smooth  and  the  starting  vector  is  calculated  by  a  method  correct 
up  to  (Xh^2),  then  the  error  in  y„  the  approximation  to  y(x,)  calculated  using  the  method  in 
Equation  2  is 


He  =  y(h) - yt  =  (-L -  f  By ^XyW) +  ^[01) +  °ihP+2)’  (53) 


provided  vectors  (3  and  y  meet  the  following  conditions: 
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1  )yTe  =  0 

2  )PTe  +  yTal  =0 


1)! 

4)^M  =  i 

;= i  P' 


ck~l 

— - +  YTcck  =0,k  =  2,...,p 


(54) 


Proof:  We  wish  to  calculate  the  term  of  0(hp+1)  in  the  expression  y(t,)  -  y,.  We  use  a 
Taylor  expansion  for  the  true  solution  about  t0  and  use  the  expression  for  the  first 
component  of  the  Nordsieck  vector  in  computing  yv  Then  we  have 


y{h)  =  Pl^p^-hl  +o(hP+2), 


We  use  the  definition  of  the  starting  vector  and  the  stage  order  condition  to  write 


y\0]=  iaijy^{t0)y  +  o(hP+1\ 

j= 0 

41]  =  y(t0+cih)  +  o(hr+l). 


Now 


^(y111)  -  /(#)  =  f(y{t 0  +  c,h)  +  0(F«))  =  f(y(t0  +  c,h))  +  0{h<’*') 

=  i(t0+cih)  +  O[h^')=  I  y!'+ll(/0)^  +  O(/.'"1). 

j=0  J- 


Then  the  term  in  h  of  order  p  +  1  in  y(tj)  -  yt  is  the  lte  and  may  be  readily  found  to  be 
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Ite  =  hp+1y(p+l\t0) 


(p  +  i)i 


7-2$ 


j= 1  1;  p! 


We  wish  now  to  find  an  estimate 


Z^+yp+i^o)  _  hpTF[Y[1]))  +  rTy[0]  +  o[hp+2). 


We  use  the  same  expressions  for  F  and  the  external  stages.  This  yields 


P  C; 


h’+1Jr*%)='ZpjZ%k 

7=1  jfc=0 


t+V‘+1)(<o)+  X  rj  I  aJi/)((0)fti  +  o(A'’-2). 

j=l  j‘=0 


The  first  relationship  comes  from  considering  terms  of  order  0.  The  second  relationship 
results  from  equating  terms  of  order  1  and  noting  that  if  some  Cj  =  0,  the  term  of  order  1  is 
produced  by  multiplying  the  term  of  order  0  in  the  Taylor  expansion  by  h  and  hence  we  use 
a  convention  where  1  appears  in  place  of  the  apparent  undefined  0°  in  the  last  formula.  The 
third  set  of  formulas  comes  from  considering  orders  2  through  p,  while  the  fourth  formula 
comes  from  setting  the  coefficients  of  terms  of  order  p  +  1  on  the  right  to  one.  ■ 

We  note  that  we  have  p  +  2  equations  in  either  case  to  determine  2p  variables.  For 
p  =  2  the  solution  is  unique,  while  for  higher  orders  there  are  additional  free  parameters. 
These  may  be  used  to  accomplish  other  purposes,  for  example  to  avoid  poles  that  might 
arise  in  the  formula. 

Error  estimation  for  variable  step  implementations  have  also  been  developed.  We 
follow  here  the  Nordsieck  approach  of  Butcher  and  Jackiewicz  (Reference  3)  but  also  note 
that  an  alternative  formulation  has  been  developed  in  a  paper  by  Jackiewicz,  Vermiglio  and 
Zennaro  (Reference  23).  We  define  hn  as  tn-tn.,,  and  8  =  hn/hn_,  and  we  seek  vectors 
|3  =  (3(8)  and  y  =  y(S)  such  that 


vVrV^Vi) = kpt(s)f(yW) + r  W-11 + o(*r2).  (55) 


It  should  be  noted  in  the  following  modification  of  a  theorem  by  Butcher  and  Jackiewicz 
(Reference  3)  that  the  error  formula  includes  the  effect  of  rescaling  and  that  error  estimation 
is  not  carried  out  by  simply  using  a  fixed  step  formula  with  a  rescaled  external  stage  vector 
as  is  done  in  interpolation.  Also  note  that  variable  stepsize  does  not  apply  before  the 
second  step.  Note  that  in  the  following  the  validity  of  die  local  stage  order  condition  is 
assumed. 


39 


NAWCWPNS  TP  8340 


Theorem  9:  Assume  that  the  method  in  Equation  2  with  p  =  q  =  r=  s,  U  =  I,  is 
implemented  in  variable  stepsize  mode  using  the  Nordsieck  technique  and  that  V  is  rank 
one,  Ve  =  e.  Then 


vTfe(i„)  =  kPt{S)f(y w) + rT(s)yl"->] + o(K*2), 


(56) 


if  p  =  j3(S)  and  y  =  y(<5)  satisfy  the  system  of  equations 


1)  yTWDV  =  0 

2) YTe  =  0 

3) (pT  +±yTWDB)e  =  0  (57) 

WT-r-^r+-brTm>BGii= o,  1  =  2,3 . P 

0-1)!  8 

5)PT^  +  -^nrTWDBGtp„  =  vT<pp. 


For  the  frequently  occurring  case  where  the  first  row  of  V  is  vT  and  the  other  rows  are  0, 
the  first  condition  simplifies  to  y  e  —  0,  eliminating  it  as  a  separate  condition. 

Proof:  We  follow  the  proof  of  Butcher  and  Jackiewicz  (Reference  3)  with  some 
modifications.  We  proceed  as  in  the  proof  of  Theorem  6.  First,  we  use  the  stage  order 
condition  and  Equation  48  to  obtain  through  Taylor  expansion, 


1). 


A‘n- 2) 
K-O'i’n-l) 
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Taylor  expansion  yields  the  relationship 


y(fn- 2) 

y(fn- 1) 

K-iy'^n-l) 

=  f 

K-iy'{tn-i) 

hPn-ly{P+%n-l). 

JCI ly{p+%n-i)_ 

and  rescaling  with  D  =  diag{ 1, \ ,  -i- , 


V 


S’S2 


1  ^ 


5  8p+] 


yields  the  relationship 


yi*n- 1) 

A*n-l) 

hn-iy'(tn-l) 

=  D 

Ky'{fn- 1) 

h^yiP+%n-i). 

hp+ly{p+%- 1). 

Substituting,  we  then  obtain 


yifn- 1) 
Ky'^n- 1) 


K+1y{p+%n-i) 


Furthermore,  using  the  relationship  in  Equation  22  and  the  rescaling  formula  in  Equation 
29,  we  obtain 


55b-1]  =  +  WDVy[n~2]. 


We  now  substitute  into  Equation  55  to  obtain 
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v>AP+VP+1)('«-i)  =  (pTG  +  yTWDBGfD) 


y{fn- 1) 

Ky'ik- 1) 

K+'yiP+1){tn-i\ 


+  yTWDVy[n~2]  +  0(hl f+2). 


Since  fixed  stepsize  D  and  Z)  are  identity  matrices,  we  may  note  that  this  expression  agrees 
with  the  comparable  expression  derived  earlier  for  fixed  stepsizes.  However,  the  term 

yTWDVy[n~1]  requires  additional  consideration.  We  note  that  the  first  column  of  G  is  0. 
Then,  examining  the  terms  created  in  multiplying  the  Nordsieck  vector  by  fiT(8)G,  we 
find  that  there  is  no  term  remaining  of  order  h°.  Similarly  since  D  is  a  diagonal  rescaling 
matrix  and  T  is  lower  triangular,  there  is  also  no  term  of  order  h°.  Then,  since 
-[n-2]  _  ey(tn_2^  +  0(hn_i)  we  must  have  yTWDVe  =  0.  For  the  case  where  the  first  row 

of  V  is  vT  and  the  other  rows  are  0  we  have  WDVe  =  e  and  the  zeroth  order  condition 
reduces  to  the  requirement  that  yT(8)e  =  0,  which  also  eliminates  the  entire  term  since 
WDV  =  WV  =  V  =  evT  for  this  case.  However,  for  more  general  V  matrices  it  has  been 
convenient  to  impose  the  sufficient  condition  yTWDV  =  0 . 

We  now  use  our  knowledge  of  the  nature  of  G  in  the  first  term  and  D  in  the  second  to 
find  that 


l fG  +  yTWDBGTD  = 


0 


(3Te  +  ±yTWDBGtl 


j~^jJTWDBGlp+i 


We  note  that  Gtl=e.  Thus  the  first  order  condition  simplifies  to 
+  ^yTWDB]je  =  0.  The  other  equations  of  Equation  37  then  follow  by  equating  terms 
of  the  same  degree.  ■ 


TECHNIQUES  FOR  OBTAINING  STARTING  VALUES  FOR  DIMSIMS 


In  general  it  is  the  derivatives  that  must  be  computed  to  obtain  starting  values. 
Techniques  will  be  illustrated  here  for  methods  of  second  order  that  may  be  extended  to 
methods  of  higher  order.  For  a  second  order  method  only  the  second  derivative  becomes  a 
problem,  since  the  first  derivative  may  be  computed  using  the  derivative  function  f  of 
Equation  1.  The  derivative  of  the  derivative  function  could  be  computed  symbolically 
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whenever  f  is  available  in  a  symbolic  form,  but  this  will  be  too  complicated  for  many 
functions  of  interest.  Otherwise,  an  approximation  for  y"  may  be  calculated,  which  needs 
only  to  be  correct  to  1st  order  to  provide  satisfactory  starting  values.  Only  one  additional 
solution  point  is  needed,  y(x0  +  h0),  where  h0  is  some  small  value  (note  that  the  first 
DIMSIM  integration  step  size  is  designated  h,).  The  Taylor  expansion  at  x0  then  yields  a 
convenient  expression  with  a  first  order  error: 


y(h ) = y(k ) + V'(?o ) + \  y\k ) + )  > 


so  we  have  the  expression 


vi-yp 

L 


yo 


+o{h o). 


(58) 


Here  we  substitute  the  approximation 


yi  =y(tl)  +  (9(/io), 


where  y,  must  be  calculated  using  a  second  order  Runge-Kutta  method  to  provide  a  second 
derivative  correct  to  O(ho).  It  may  be  noticed  that  an  approximation  for  y(tj)  is  obtained  as 
part  of  the  starting  procedure.  However  it  is  the  external  stages  that  are  needed  and  they 
will  be  calculated  at  t,  using  the  DIMSIM.  New  function  evaluations  will  be  required,  even 
if  the  value  of  c,  is  0  and  h0  is  actually  used  as  the  initial  step  size. 

Although  adequate  starting  values  are  obtained  with  a  first  order  estimate  of  the  second 
derivate,  as  was  shown  in  Theorem  8,  in  order  to  obtain  a  reliable  error  estimate  for  the 

first  step  for  a  second  order  DIMSIM,  it  is  necessary  to  obtain  an  o[hq )  approximation  for 

y"(t0) .  This  may  be  done  using  a  3rd  order  Runge-Kutta  method  and  either  two  calculated 

points  or  one  calculated  point  and  a  functional  evaluation,  as  follows.  Assume  we  calculate 
two  points  yj  and  y2  from  (t0,y0)  using  a  3rd  order  (explicit)  method.  Note  that  if  both  are 
calculated  from  the  starting  point  there  will  be  no  stability  problems,  no  subsequent  steps 
are  taken  with  the  method.  We  can  express  both  in  terms  of  Taylor  expansions  of  y(t0+hj) 

4  4 

and  y(x0+h2)  to  0(h  ),  since  these  are  the  accurate  with  error  0(h  ).  Let  h]  =  h{/2  and 
h2  =  h0.  Then  we  have: 
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y>=yo+bi+iy'o+iyS'+0(h°) 
=  y0  +  «  +4>o+f  ?o'+0(*o4). 


We  can  eliminate  the  third  derivative  term  between  these  two  equations  and  come  up  with 
the  following  expression  for  the  second  derivative: 


y'o  =  fe  -  8ti  +  7y0 + ihoy'o ) + o{$  )• 


(59) 


Alternatively  we  can  calculate  a  yj  at  t0+h0,  do  a  function  evaluation  there  to  determine  its 
first  derivative  y{,  and  use  Taylor  expansions  to  develop  a  different  formula,  also  correct  to 
second  order: 


yi=y0 + W  +  4  >? + 4  /o+ o(hS ) 
y!  =  yo+W+4yS ’+o(4\ 


Multiplying  through  the  second  equation  by  h0  we  have  an  error  term  that  is  4th  order,  the 
third  derivative  can  again  be  eliminated  and  we  arrive  at  the  formula: 


y'o  =  7%(3(vi  -  To)  -  ho(y{  +  2 To))  +  o(f$ ) . 


(60) 


Limited  testing  seemed  to  indicate  that  Equation  59  produced  somewhat  better  results. 
However,  the  second  method  requires  fewer  function  evaluations.  But  a  great  advantage  to 
either  approach  is  that  they  yield  a  convenient  estimate  for  the  third  derivative,  and  this  may 
be  used  to  provide  a  very  accurate  a  priori  error  estimate  enabling  optimal  step  size 
selection.  This  will  be  explained  in  more  detail  in  the  next  section. 


STEP  SIZE  SELECTION  STRATEGY 

The  adaptive  approach  employed  for  each  step  is  to  use  the  error  estimation  techniques 
outlined  above  to  obtain  an  estimate  of  the  error  generated  with  the  step.  If  the  error 
exceeds  the  tolerance  the  step  size  is  halved,  the  external  stage  vector  is  rescaled  using 
Equation  29,  and  the  step  is  repeated.  This  continues  until  a  result  is  obtained  within 
tolerance  or  the  maximum  allowable  number  of  attempts  is  exceeded,  which  terminates  the 
integration  with  an  error  message.  On  the  other  hand,  if  a  step  is  successful,  a  new  step 
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size  is  calculated  for  the  next  step.  This  is  done  with  the  standard  formula  provided,  for 
example,  by  Hairer,  Norsett,  and  Wanner  (Reference  18).  The  optimal  step  size,  hopt,  will 
produce  an  error  equal  to  the  tolerance.  The  error  is  assumed  for  both  the  previous  step 
and  the  next  step  to  be  of  the  form  of  the  left  hand  side  of  Equation  55  with  approximately 
the  same  derivative  factor.  Then  we  estimate 


^opt  ^ 


1  V+1 

err  > 


There  are  some  issues  that  must  be  considered  in  changing  step  size,  even  after  a 
computation  of  optimal  step  size  has  been  carried  out.  First,  the  amount  of  work  that  is 
necessary  to  complete  an  integration  is  minimized  if  the  longest  possible  step  sizes  are 
used.  On  the  other  hand,  failed  steps  are  expensive  since  all  the  step  calculations,  including 
rescaling  external  stages,  derivative  evaluations,  and  error  estimation  must  be  repeated. 
Therefore,  since  error  estimates  are  simply  estimates,  a  safety  factor  is  always  utilized, 
usually  chosen  to  be  0.8  or  0.9.  Furthermore,  large  changes  in  optimal  step-size 
calculations  are  an  indication  of  rapid  changes  in  the  solution  to  the  system  being  solved 
which  lessens  the  value  of  error  estimates,  and  so  a  maximum  step-size  increase  ratio  is  set. 
Finally,  error  estimates  and  integration  behavior  deteriorate  with  frequent  step-size 
changes,  so  it  is  desirable  once  a  step  size  has  increased  to  prohibit  further  increases  for  a 
few  steps.  There  are  then  three  parameters  to  be  determined  heuristically  that  can 
significantly  affect  the  performance  of  a  solver,  the  safety  factor,  the  maximum  step-size 
increase  ratio,  and  the  number  of  successive  steps  after  a  step-size  change  in  which  step 
size  is  kept  from  increasing. 

Integration  at  the  right  end  point  may  proceed  in  a  few  different  ways.  If  the  interval 
from  the  last  mesh  point  to  the  end  point  is  small,  extrapolation  may  be  used. 
Alternatively,  the  solver  may  integrate  past  the  end  point  using  the  step  length  calculated 
from  the  error  estimate  and  interpolate  back.  Finally,  the  right  end  point  may  be  chosen  as 
the  final  mesh  point  for  the  last  integration  step.  This  third  alternative  is  used  here  because 
it  is  most  appropriate  for  waveform  relaxation.  Extrapolation  does  not  yield  an  interpolant 
for  the  final  interval,  which  is  the  actual  required  output  with  waveform  relaxation.  Also, 
in  waveform  relaxation  the  derivative  function  is  undefined  past  the  end  of  the  current 
window  since  interpolants  for  variables  associated  with  other  subsystems  of  the  overall 
problem  have  not  yet  been  determined.  Thus  the  third  alternative  is  the  technique  chosen 
for  this  work. 

Choosing  a  suitable  size  for  the  first  step  has  been  almost  as  much  of  an  art  as  a  science 
and  approaches  are  typically  taken  based,  to  a  considerable  extent,  on  heuristics.  It  is 
desirable  for  accuracy  to  choose  a  very  small  first  step  and  then  use  the  error  estimator  to 
determine  step-size  changes  for  subsequent  steps.  However,  too  small  a  first  step  will 
result  in  too  many  small  steps  as  the  integration  progresses  while  the  step  size  is  increasing. 
Three  different  approaches  were  examined  in  developing  DIMSIM  solvers.  Shampine  and 
Gordon  (Reference  24)  used  the  following  selection  for  initial  step  size: 
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h  =  min 


1 

’4 


0.5*toZ 

|/('o>>’o)i 


1 

2 


(61) 


Here  the  input  h0  value  is  a  user-supplied  estimate  to  prevent  too  large  a  value  from  being 
used,  and  tol  is  the  specified  tolerance.  The  choice  comes  from  estimating  the  error  of  a 
first  order  method  as  h  times  the  error  of  a  zero  order  method,  and  then  calculating  a  step 
size  to  produce  an  error  of  half  the  tolerance.  This  step  size  is  then  divided  by  4  to  produce 
a  conservative  value.  Their  approach  is  used  in  a  linear  multistep  solver  in  which  the  first 
step  is  of  first  order,  which  renders  it  less  suitable  for  a  higher  order  DIMSIM  solver. 

A  more  sophisticated  approach  suitable  for  higher  order  solvers  is  presented  by  Hairer, 
Norsett  and  Wanner  (Reference  18)  and  originally  developed  by  Gladwell,  Shampine,  and 
Brankin  (Reference  25).  The  local  truncation  error  is  assumed  to  be  of  the  form 


he  *  Chp+ly^p+l\tQ). 


(62) 


They  then  recommend  the  following  process.  The  preliminary  value  is  the  step  size  that 
produces  an  Euler  step  yielding  a  solution  change  of  1%.  This  is  used  to  then  obtain  an 
estimate  for  the  norm  of  the  second  derivative.  The  larger  of  the  first  and  second 
derivatives  is  used  as  an  estimate  for  Cy(p+1)(t0),  and  a  value  of  h  is  chosen  to  produce  an  lte 
of  around  1%.  Various  threshholds  are  set  to  avoid  bad  choices  for  more  exceptional 
cases.  In  summary,  these  steps  are  followed: 


1.  Let  dQ  =  ||yo||  and  dx  =||/(Wo)||’  where  Nl  =  , 

scale  factors,  here  sc—Atolj+lyJRtoli  with  Atol  and  Rtol  vectors  of  absolute  and  relative 
error  tolerances. 

2.  If  d0  or  d;  is  less  than  10'5,  set  h0=  10'6,  otherwise  let  h0=  0.01  (d,/d,). 

3.  Let  y,  =  y0+  h0f(t0,y0)  and  find  f(t0+  h^). 

4.  Let  d2  =i||/(/0  +^,yi)-/('o>}’o)||- 

5.  If  max(Jj,J2)^  10-15  then  set  hx  =  max^lO-6,/^  -10~3j,  otherwise, 

0.01 

max^j,  J2) 


,  n( 

ly 

n  ", 

i=lV 


Zj_ 

SC; 


and  sc  is  a  vector  of 
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6.  Use  starting  step  size  h  =  min(100h0,h1). 

This  was  the  first  method  that  was  tried  for  choosing  the  initial  step  size. 

An  alternative  approach  with  considerable  promise  has  been  devised  as  part  of  this 
research  effort.  For  a  5th  order  DIMSEM,  for  example,  we  would  proceed  as  follows:  Any 
small  step  may  be  used  with  an  explicit  6th  order  Runge-Kutta  integrator  to  produce  3 
values  for  the  solution  at  t0  +  h0,  t0  +  2h0,  and  t0  +  4h0.  Using  two  of  the  three  possible 
derivative  function  evaluations  to  produce  values  for  the  derivative,  there  are  5  equations  in 

5  unknowns,  producing  values  for  the  /zoy^(t0)  for  k  ranging  from  2  to  5  as  needed  for 

the  starting  vector,  plus  a  value  for  ^o;/6^o)>  which  may  be  used  with  the  existing 
formula  to  provide  an  accurate  a  priori  error  estimate  for  the  first  step  before  the  DIMSIM  is 
even  applied,  using  Equation  53.  This  may  be  used  with  the  tolerance  value  to  provide  a 
simple  formula  for  the  appropriate  step  size  for  the  first  step.  A  safety  factor  of  perhaps  2 
should  be  provided  based  on  heuristics  from  further  tests.  The  starting  vector  calculated 
using  h0  would  only  need  rescaling,  prior  to  actual  use.  This  will  avoid  the  usual  step  size 
buildup  at  the  beginning  and  improve  accuracy.  Furthermore,  one  might  utilize  alternative 
error  formulas  for  solvers  of  different  orders,  with  no  additional  derivative  function 
evaluations  required,  to  obtain  an  optimal  order  selection  for  the  problem,  at  least  in  the 
region  of  the  first  step. 

Some  problems  encountered  with  early  testing  of  this  approach  indicate  that  it  is 
important  to  utilize  a  suitable  h0  value.  It  should  be  noted  that  the  most  difficult  value  to 

calculate  in  the  Nordsieck  vector  will  be  the  term  of  0(h£+l)  if  h0  is  small,  because  it  is 
calculated  as  a  linear  combination  of  terms  that  may  be  near  1.  Thus  machine  precision  and 
roundoff  error  may  become  significant.  The  problem  of  machine  precision  is  greatly 
simplified  when  a  quadruple  precision  data  type  (real*  16)  is  available.  Furthermore,  if  h0 
is  too  large,  the  Runge-Kutta  integration  will  not  be  sufficiently  accurate.  The  estimate  for 
a  suitable  value  of  h0  may  be  obtained  from  the  heuristic  relationship 


1 0  •  eps  =  scfac  ■  h£+l , 


where  eps  is  the  machine  epsilon  and  scfac  is  a  scale  factor  reflecting  the  scaling  of  the 
problem.  A  reasonable  choice  of  scfac  might  be  1  or  the  minimum  of  y0  and  y0’  with  the 
maximum  taken  if  one  of  the  two  is  very  small  (an  approximation  to  0).  This  assumes  that 
the  (p  +  l)st  derivative  is  around  the  same  size  as  scfac.  The  idea  is  to  choose  hQ  so  that  the 
smallest  and  most  difficult  term  to  calculate  is  on  the  order  of  ten  times  the  machine  epsilon. 
If  h0  is  smaller,  roundoff  error  will  create  problems,  while  if  h0  is  too  large  the 
approximation  accuracy  will  suffer.  This  process  is  more  extensively  illustrated  as  it  is 
used  for  obtaining  starting  procedures  in  the  Techniques  for  Obtaining  Starting  Values  for 
DIMSIMs  section. 
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ALTERNATIVE  REPRESENTATIONS 

The  DIMSIM  solution  equation  and  the  rescaling  formula  provide  relationships  among 
current  external  stage  vectors,  previous  external  stage  vectors,  and  derivative  vectors  at 
current  or  previous  internal  stage  points.  These  enable  alternative  representations  of 
formulas  for  such  things  as  error  estimates,  intepolants,  Nordsieck  vectors,  and  rescaling. 
These  alternative  formulas  are  mathematically  completely  equivalent  but  may  have 
somewhat  different  roundoff  properties.  In  some  cases  significant  differences  for 
waveform  relaxation  implementations  in  memory  requirements  and  in  message  passing 
volume  for  parallel  computing  will  result  from  the  choice  of  representation,  as  will  be 
described  below.  Alternative  representations  are  also  helpful  in  shortening  recursions  in 
deriving  stability  regions  for  predictor-corrector  implementations,  and  this  will  be  used  in 
the  report  to  follow  on  implicit  DIMSIMs.  The  following  relationships  are  not  exhaustive 
but  are  provided  as  examples  of  useful  forms. 

We  again  write  the  second  equation  of  Equation  2,  using  an  alternative  expression  for 
V: 


yM  =  fcBFjyM)  +  evTy^n~^. 


(63) 


We  can  then  write,  for  nonsingular  B,  (implicitly  assumed  in  this  section  wherever  B'1 
appears) 


-  B~levTy^n~l\  (64) 

Alternatively,  we  may  use  the  rescaling  formula  in  Equation  29  to  express 
=  hn^WDBF[Y[n~l]^  +  WDVy[n~2]. 

We  can  then  apply  these  as  follows: 

a.  Error  Estimates:  The  fixed  step  formula  has  been  expressed  in  Equation  50  as 
vTle(tn)  =  hpTF^n^  +  yTy[n~l]  +  0[hp+ 2) , 
and  we  may  produce  alternative  representations  as  described  in  the  following. 
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Corollary  to  Theorem  7:  If  the  fixed  step  method  in  Equation  2  has  stage  order  q  equal 
to  the  order  p  and  V  is  a  rank  one  matrix  such  that  Ve  =  e,  then  vTle(xn)  can  be  estimated  by 
Equation  50  and  either  of  the  formulas 

vTle(tn )  =  +  0[hp+iy  (65) 


and 


vTle(tn)  =  Yi  +  rh[n  1]  +  0{hp+2), 


(66) 


rp  rp 

where  p  and  y  are  given  by  the  formulas  of  Equation  51,  p,=  p,  y32  =  /  5,  and  y,  and  y2 
are  given  by 

rf =Ptb~\ 

(67) 

yT2=yT-f}rB-W. 


Proof:  To  obtain  Equation  66  we  substitute  Equation  64  into  Equation  50  to  obtain 


vTle(tn)  =  (i r(s-V"l  -  B-Win~x  1)  +  rTyln~']  +  0(hP*2), 

and  the  result  follows  immediately.  Equation  65  results  from  substituting  Equation  63  for 
the  previous  step  into  Equation  50  to  obtain 

vTle(tn)  =  h{3TF{Y[n]])  +  YT{hBF[Y[n^  +  evTy^n~2^. 


The  extra  term  yTevTy^n  is  0  since  yTe=0  from  Equation  51.  ■ 

We  have  similar  relationships  for  variable  step  size.  Note  that  5  refers  to  8n  in  the 
following. 

Corollary  to  Theorem  9:  If  the  variable  step  method  in  Equation  30  has  stage  order  q 
equal  to  the  order  p  and  V  is  a  rank  one  matrix  such  that  Ve  =  e,  then  vTle(xn)  can  be 
estimated  by  Equation  56  and  any  of  the  formulas 
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vrM'»)=*»Ar  (A)^^“1)+A,-ifer(A.A-i)f(^""11)+o(^+2),  (68) 


and 


/%) = rl  (<%w + jf  WM1 + o^2). 


(69) 


where  p  and  y  are  given  by  the  formulas  of  Equation  57,  P,=  P,  fa  ~  7 ^ (<5rt )  WD(8n_i ) B , 
with  V  in  standard  form,  and  y,  and  y2  are  given  by 

rt  =  PT(5)B-\ 

(70) 

yT2=yT(8)-J3T(8)B-W. 

Proof:  To  obtain  Equation  69  we  substitute  Equation  64  into  Equation  56  to  obtain 

vrfe(t„) = hf  (S)(fi-y”)  -  CTV“1]) + yT(S)in~'], 

and  the  result  follows  from  identification  of  coefficients.  Equation  68  results  from 
substituting  Equation  29  into  Equation  56  to  obtain 

vTle(tn)  =  hnf(Sn)F[Y^YrT{Sn{wD{Sa_l)^n-l^  +  WD^"-^y 


Assuming  the  standard  form  of  V  with  a  first  row  of  vT  and  all  the  other  rows  identically  0, 
we  may  note  that  WDV  =  V  =  evT .  This  results  in  the  extra  term  yT (8)evT y^n~2\  but 
yT(8)e  =  0  from  Equation  57.  ■ 

A  similar  approach  can  be  used  for  the  first  step  error  estimate. 

b.  Interpolants: 

The  most  payoff  from  using  alternative  formulas  can  be  expected  here,  since 
interpolants  over  some  number  of  steps  are  stored  and  passed  between  parallel  processors 
as  they  are  used  to  represent  solutions  over  a  window  in  waveform  relaxation.  The 
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formula  derived  above  requires  storing  and  passing  both  external  stage  vectors  for  each 
step  within  a  window  and  the  vectors  of  derivatives  at  internal  stage  points,  at  each  grid 
point  requiring  storage  and  passing  of  2Mp  numbers,  where  M  is  the  number  of  equations 
in  a  subsystem  and  p  is  the  order  of  the  DIMSIM.  The  requirement  for  a  Nordsieck 
interpolant  is  only  M(p  +  1),  but  in  a  representation  using  only  external  stage  vectors  this 
number  becomes  just  Mp.  We  have  the  following  alternative  formulations.  It  should  be 
noted  carefully  which  external  stage  vectors  are  rescaled. 

Corollary  to  Theorem  5:  If  a  DIMSIM  has  order  p  and  stage  order  q  =  porq  =  p-  1, 
then  a  Butcher-Jackiewicz-type  interpolant  77  given  by  either  of  two  equivalent  forms 
approximates  y  with  uniform  order  p  and  is  continuous  on  the  whole  interval  of  integration 
if  and  only  if  conditions  of  Equations  35  and  36  are  met  and  the  interpolant  is  given  by 
Equation  34  or 


n{tn- 1  +dhn)  =  YoMy[n] +7o,2(0)5>["_11> 


(71) 


where 


7o,i(e)  =  Po(d)B~1’ 


(72) 


7o,2(0)  =  7o(0)-ft)(0)£  XwT ■ 

Proof:  This  follows  from  identification  of  coefficients  after  direct  substitution  of  Equation 
64  into  Equation  34  to  obtain 

+  6K)  =  A )(0)(^_1y^  -  B~levTy^n~1^  +  7o(0)^n_1]-  ■ 


For  the  Nordsieck  vector  we  have  the  following. 

Corollary  to  Theorem  3:  If  the  method  in  Equation  2  satisfies  Equation  24,  then  the 
Nordsieck  vector  y  at  tn  may  be  approximated  to  0(hp+1)  using  either  the  equivalent 
formula  of  Equation  22  or 

5y„)=K>w+'y[""1,>  p3) 


where 
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Vx  =  BB-\ 

(74) 

V2  =  V-BB~1evT. 


Proof:  This  follows  from  identification  of  coefficients  after  direct  substitution  of  Equation 
64  into  Equation  22  to  obtain 


y(t„)  =  b(b~ V"1  -  B~levTin~l]^  +  Vy[n_1] .  ■ 


c.  Rescaling:  Since  formulas  used  above  include  some  rescaled  and  some  unrescaled 
external  stage  vectors,  the  ability  to  rescale  using  only  external  stage  vectors  is  important. 
We  note  that  if  we  can  rescale  using  Equation  29,  we  can  also  rescale  using 


y[n~1]  =  WDnBB~Xy[n~X]  +  (/  -  WDnBB~X  )evTy[n~2] .  (75) 


This  follows  immediately  from  direct  substitution  of  Equation  64  into  Equation  29  to  obtain 


=  Vt  WD{8n)B[B~ly^n~^  -  B~levT +  WD{8n)elvT y[n  2\ 


where  a  standard  form  for  V  is  assumed  and  I  represents  the  identity  matrix  of  appropriate 
dimensionality. 


FIRST  APPROXIMATELY  SAME  AS  LAST  (FASAL)  MODE 

Runge-Kutta  methods  for  which  the  first  stage  of  the  new  step  equals  the  last  stage  of 
the  previous  step  have  been  know  for  some  time  (see,  for  example.  Reference  18),  and  the 
property  is  called  First  Same  As  Last  (FSAL).  Its  primary  use  is  for  error  estimation,  as 
for  example  with  Dormand-Prince  pairs  (Reference  26),  where  an  embedded  method  of 
higher  order  is  created  by  adding  an  extra  stage  with  A  coefficients  identical  to  the  b 
coefficients  of  the  first  method.  The  higher  order  result  is  actually  used  to  continue  the 
integration  in  this  case  but  seven  stages  are  utilized  for  a  5th  order  method.  No  additional 
function  evaluations  are  needed  for  successful  steps  and  the  work  is  the  same  as  though 
only  6  stages  were  used  because  the  final  stage  is  identical  to  the  first  stage  of  the  following 
step  (FSAL)  which  would  have  to  be  evaluated  anyway.  Because  of  the  high  stage  order, 
all  types  of  DIMSIMs  for  which  ct  =  0  and  cp  =  1  and  the  local  stage  order  condition 
(Equation  41)  is  true  have  a  First  Approximately  Same  As  Last  (FASAL)  property,  and  this 
can  be  used  to  save  at  least  one  internal  stage  evaluation  on  every  step  after  the  first.  As 
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discussed  above  in  conjunction  with  error  estimation,  theoretical  work  concerning  the  effect 
of  global  error  on  local  stage  error  remains  to  be  done,  but  experimental  evidence  for  the 
validity  of  this  assumption  is  abundant. 

The  local  stage  order  condition  requires  that 


Yp*  ^  “  y({n-2  +  +  °[K-1  )’ 


where  here  y  denotes  the  local  solution  of  the  ODE  through  the  indicated  initial  point,  while 


4n]  =  +  o(h%+1), 


where  y  again  is  a  local  solution,  and  hD  =  hn_a  and  tn.,  =  tn_2  +  hn_r  Note  that  the  order 
condition  implies  that 


yi^n— l)  —  yifn-2  2’^n— 1  )’ 


and  thus  we  have,  with  h  =  hn=  hn.1; 


y[n]  _  y[n- 1]  +  o^hp+ 1  j 


Since  what  is  used  is  the  always  the  derivative  multiplied  times  the  step  size,  and 


+ o[h”*'))  --  v(tn-i,rin]) + o(hp*2). 


it  is  possible  to  use  f{y^  in  place  of  carrying  out  the  function  evaluation  for  explicit 
methods,  or  even  the  nonlinear  equation  solving  required  in  implicit  methods,  for 


It  should  be  noted  that  higher  order  terms  (0(hp+2))  will  be  changed.  In  some  cases 
they  will  be  fortuitously  decreased,  but  in  others  they  will  be  increased.  The  smaller  the 
error  constant  the  greater  the  impact  that  these  terms  will  have.  But  in  general  the  error 
estimation,  rescaling,  and  interpolation  for  the  original  DIMSIM  should  be  essentially 
unchanged. 
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But  the  most  significant  impact  will  be  on  stability.  It  should  be  noted  that  when  this 
implementation  is  used,  the  pth  internal  stage  from  the  previous  step  is  carried  over  into  the 
next  step  along  with  the  external  stage  vector,  and  thus  in  a  sense  this  becomes  a  “two-step 
DIMSIM.”  The  modified  method  may  be  written  in  standard  General  Linear  Method 

AAA  A 

(GLM)  notation  by  producing  A,  B,  U ,  and  V  matrices  from  the  A,  B,  U(=  I),  and  V 
from  the  original  DIMSIM,  and  by  enlarging  the  external  stage  vector  length  by  1  to 
r  =  p  +  1.  This  becomes  a  GLM  with  p  =  q=s,  r  =  p+l.  The  first  component  of  the 
external  stage  vector  becomes  the  pth  internal  stage  that  is  carried  along,  while  the  others 

are  unchanged.  We  set  A  =  A  except  for  the  first  row  which  becomes  0;  this  is  unchanged 

A, 

for  explicit  methods.  U  is  an  augmentation  of  I  produced  by  adding  a  0  second  column, 
and  hence  is  px(p  +  1).  B  is  produced  by  adding  to  B  the  last  row  of  A  as  a  new  first  row 

A 

and  is  (p  +  l)xp.  Finally,  V  is  produced  by  adding  to  V  a  new  first  row  and  first  column, 
with  1  for  the  last  element  of  the  first  row  and  Os  elsewhere  for  these  new  elements.  A 
tableau  then  appears  as 

0  10  0-0' 

A2  0  0  1  : 

:  :  :  \  0 

Ap  0  0  •••  0  1  ’ 

Ap  0  -  1 

B  0  F 

where  Ak  is  used  to  denote  the  kth  row  of  matrix  A. 

The  original  V  becomes  the  lower  right  submatrix  and  the  new  matrix  is  (p  +  1)  x  (p  + 
1).  Then  we  may  write  for  a  single  ODE, 


A  U 

A  A 

B  V 


yW  =  mf(fw)  +  Uytn~l\ 
yW  =  +  Vy[n-1l 


(76) 


(It  should  be  noted  here  that  the  external  stages  are  redefined,  although  notation  has  not 
been  changed,  and  that  actual  implementation  will  not  involve  separate  evaluations  or 
function  evaluations  for  the  repeated  stages.)  Then  the  new  stability  matrix  becomes, 
through  application  to  the  test  problem  y’  =  \y  and  setting  h?i  =  z, 


M- 


■  zB{l  -  zA  j 


-1 


u+v. 


(77) 
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This  is  now  (p  +  1)  x  (p  +  1)  and  in  general  does  not  produce  the  same  stability  region  as 
the  original  DIMSIM.  Thus  if  FASAL  implementation  is  desired,  a  method  should  be 
sought  that  optimizes  the  FASAL  stability  region  rather  than  that  of  the  DIMSIM,  and  in 
general  it  should  not  be  expected  that  FASAL  implementation  for  a  given  DIMSIM  should 
produce  a  favorable  region.  This  may  be  used  in  various  ways,  including  development  of 
new  methods  especially  for  FASAL  implementation,  or  in  variable  implementation  where 
FASAL  is  used  until  it  runs  into  stability  problems.  Typical  stability  region  plots  will  be 
shown  below  as  the  use  of  specific  methods  are  discussed. 


DEVELOPING  EXPLICIT  DIMSIM  ODE  SOLVERS 


Explicit  solvers  have  been  developed  based  on  Type  1  DIMSIMs  derived  by  Butcher 
(Reference  1)  and  by  Butcher  and  Jackiewicz  (References  6,  7,  and  8).  These  methods  are 
designed  to  have  the  same  desirable  stability  regions  as  Runge-Kutta  methods  but  break 
Runge-Kutta  order  barriers  at  orders  above  4  and  stage  order  barriers,  thus  reducing  the 
number  of  function  evaluations  required  per  integration  step.  These  DIMSIMs  are 
designed  to  have  a  number  of  stages  (the  number  of  function  evaluations)  equal  to  the  order 
of  the  method,  while  explicit  Runge-Kutta  methods  require  one  extra  stage  for  orders  5  and 
6,  two  extra  for  order  7,  and  at  least  three  extra  for  orders  8  and  higher.  In  fact,  order  10  is 
the  highest-order  explicitly  constructed  explicit  method  found  so  far,  and  the  fewest  stages 
required  of  any  10th  order  method  is  17  according  to  the  current  (1993)  revision  of  Hairer, 
Norsett,  and  Wanner  (Reference  18).  Thus  the  recent  announcement  of  the  discovery  of  an 
8th  order  Type  1  DIMSIM  using  8  stages  and  with  good  stability  properties  by  Butcher, 
Jackiewicz  and  Mittelmann  (Reference  6)  is  of  great  significance. 

The  development  of  a  solver  involves  derivation  and  testing  of  additional 
implementation  parameters  to  provide  for  rescaling  at  step-size  changes,  error  estimation, 
interpolation,  and  starting  procedure.  A  solver  also  requires  software  design  and 
implementation.  In  this  case  the  final  codes  were  written  in  FORTRAN  77. 


IMPLEMENTION  PARAMETERS  FOR  A  SECOND  ORDER  EXPLICIT 
DIMSIM 

Although  second  order  does  not  provide  sufficient  accuracy  to  be  broadly  useful,  the 
simplicity  of  the  small  number  of  parameters  enables  development  and  convenient 
illustration  of  techniques  applicable  to  higher  orders. 

Butcher’s  second  order  Type  1  method  in  Equation  15  was  utilized  to  develop  a 
DIMSIM  variable  step-size  solver.  For  convenience  the  tableau  is  reproduced  here. 


"A  TJ~ 
B  V 


'00  1  O' 

2  0  0  1 

5.  1  11’ 

4  4  2  2 

1-111 
4  4  2  2. 
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and  c  =  [0,1  ]T.  Butcher  (Reference  1)  also  found  matrices  W,  B,  and  V  (see  above  in 
the  Introduction  to  DIMSIMs  section  for  the  relevant  definitions).  These  are 


W  = 


1  0  O' 

1  i. 


i 

4 


B  = 


0 

-1 


1 

2 


1 

2 


v= 


0  0 
0  0 


(78) 


The  latter  two  were  determined  using  a  free  parameter  g  to  obtain  desirable  step-size 
changing  zero  stability,  determined  from  the  eigenvalues  of  the  matrix 


WDV  = 


1 

2 

\  -  Sg  +  S2g 


1 

2 


"1 

r 

2 

2 

1 

2 

1 

2 

=  v. 


for  g  =  0,  which  has  eigenvalues  0  and  1,  independent  of  8.  Note  that  standard  form  is 

used  for  V  in  this  case.  The  stability  polynomial  M(z)  =  1  +  z  +  z2/2,  is  the  same  as 
for  a  two-stage,  second-order  explicit  Runge-Kutta  method,  and  thus  this  method  has  the 
same  well-known  stability  region,  including  the  interval  [-2,  0]  along  the  real  axis. 

Butcher  (Reference  1)  also  derived  estimates  for  the  local  truncation  error  (after  the  first 
step)  as 


"e  -  rfs(ihf^  -  - +  tH""11 


(79) 


which  reduces  for  fixed  step  (5  =  1)  to 
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lie  =  lhf{Y2)  -  \hf(Y,) + j(yj”-1]  -  4""11)- 


For  the  first  step  error  estimate  we  apply  Theorem  8,  Equation  53  and  calculate 


1  p  ~  c-j  . 

-z.vi-i- 


(P  +  1)!  >1  V  P! 


The  simultaneous  equations  to  satisfy  from  Equation  54  are 

71+72=° 

A+/*2-72=° 

&+i72  =  ° 

2^2  =1 


Then  pT  =  [-6,2],  yr  =  [4,-4],  and  the  first  step  error  estimate  becomes 

= &  ¥& )  -  i  ¥(? ) + i  (?’  -A01)- 


(80) 


This  is  exactly  one-half  of  the  fixed-step  error  estimate  used  for  subsequent  steps. 

The  final  computation  involves  the  derivation  of  an  interpolant  using  Theorem  5.  We 
need  to  find  the  12  coefficients  of 


ftfa)  = 


Ano  +  Ani0  + A)i2#2 

ft 020  +  A)21^  +  fto22^2 


and 


7oTW  = 


7oio  +  7on^  +  7oi2^2 
.7020  +  7o2i^  +  7  022^ 2 
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where  6  e  [0, 1]  is  the  interpolation  parameter.  We  find  immediately  from  the  continuity 
condition  |30(O)  =  0  that  we  must  have  (3010  =  P020  =  0.  The  other  continuity  conditions 
11(3.3),  when  combined  with  the  compatibility  conditions  (see  Implementing  DIMSIMs 
section),  yield  the  equations 


yQ(0)B=  f/oio+f/020  4/010  ~  4/020]-  A>(1)-[f  i]> 


which  yields  Y01o=l5  Yo2o=0>  and 


/o(0)V,  =  ._i(/oiO  +  /02o)  i(/010+/020)]“  /oC1)-  2  2]’ 


which  is  consistent  with  the  previous  result  but  yields  no  additional  information.  The  basic 
compatibility  conditions  then  tell  us  that 


f  _  A)ll+A)12 
|_A)21  +  A322J 

2  _  /010  +/011  +/012 

.iJ_L/020  +/021  +/022j 

or,  using  the  values  just  calculated,  the  second  vector  equation  becomes 
~j  _  /011+/012 

.1  J  |/021  +/022.‘ 

The  order  condition  now  becomes 
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We  expand  the  exponentials  about  z  =  0  and  equate  coefficients  of  powers  of  z.  This  gives 
us  the  following  three  equations: 


0(1) :  1  +  Ton#  +  7oi2^2  +  7021^  +  7o22^2  - 1 


0{z)  :  Pq\\0  +  A)1202  +  Po21@  +  Po22&2  ~  Yo21&  ~  Y 022^  ~  ® 
°(z2)  '■  A)21^  +  P022O2  +  2  7021^  +  J&2  =  —■ 


This  gives  us  7  equations  for  8  variables,  but  we  clearly  must  seek  to  eliminate  the 
parameter  dependence  on  0.  We  use  the  compatibility  conditions  to  eliminate  yon,  y021, 
(3on,  and  p021.  Our  0(1)  condition  then  reduces  to 


-(1  -  0)/oi2  -  0  -  0)7022  - 


which  is  used  to  eliminate  y022  when  the  condition  0  *  1  is  considered.  The  0(z)  condition 
then  yields  the  equation 


i-0-  2^22(1  -0)-  7022(1  -0)~o. 


This  similarly  enables  elimnation  of  y022.  Finally,  the  0(z2)  condition  is  applied.  It 
reduces  to 


1-O-/3o12(1-O)-3/3o22(1-0)  =  O. 


We  then  can  eliminate  P012,  leaving  us  with  the  following  one-parameter  family  of 
coefficients: 
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7oio  - 

Toil  =  2  _2A>22> 

7012  =  -1  +  2/?02’ 

7020  =  0, 

7021  =  _i  +  2^022’ 

7022  =  1  -  2^022  > 

A)10  =  0’ 

A)ll  =  7  +  2A)22  > 

A)12  =  l-3^022, 

A)20  =  0> 

2*021  =  i~  A)22> 

Po22  ~  P022  • 

We  let  (3022=  0  to  obtain  the  Butcher-Jackiewicz-type  interpolant 

n(*n- 1  +  eft) = ft(ie + e2)/(iiw) + ftf /(if1) 

+  (^l+f.e^»-i]  +  ^£  +  e^M] 


Some  possible  starting  procedures  were  derived  above  in  the  chapter  Implementing 
DIMSIMs  as  an  illustration  of  techniques  generally  applicable  to  DIMSIMs.  For  our 
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explicit  second  order  solver  we  use  an  order  3  explicit  Runge-Kutta  solver  to  obtain  an 
approximation  yl  to  y(to+h0)  accurate  to  O(h04).  Then  we  have: 


yi=y0  +  h0y'0+^ryo  + 


& 

6 


Yi'+o(h o4). 


y;=y’o+^yS+^ 


yo ’+o(h$). 


Here  y[  =  f(t0  +hQ,yl).  Eliminating  variables  we  solve  these  two  simultaneous  equations 
and  obtain 


yS=^yi-yo)-U-y{)+o(hi), 

y'o  =  tt(2(>'o  -  ^t) + y'o + yi) + )• 

"0 


Note  that  the  third  derivative  is  not  needed  for  the  Nordsieck  vector  but  is  useful  in 
selecting  initial  step  size.  Then  the  Nordsieck  vector  at  ^  is  given  by 


yM = 

o) 

— 

VtWo) 

avm. 

6{yl-y0)-4y'0-2yi_ 

Note  that  this  depends  on  the  step  size.  Once  a  correct  step  size  hj  is  chosen,  this  vector  is 
rescaled  using  the  matrix 


D  = 


1  0 
0  8 
0  0 


where  8  is  the  ratio  — .  Multiplication  by  W  then  yields  the  appropriate  starting  vector. 

K 

The  error  in  the  first  step  is  given  by  Theorem  8  to  be 
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Ite  i  = 


'  ^  Cp 

V  '=1  Z'  J 


y{3W  +  °{hf)  =  i^3(2(>o  -  yi) + h  {y'o + tf))  +  o(p) . 


If  we  set  this  error  to  half  the  tolerance  t  and  use  a  norm  to  include  the  possibility  of 
systems  of  equations,  we  then  may  write  a  conservative  but  hopefully  accurate  choice  for 
initial  step  size  to  be  8h0,  where 


8  = 


2t 


2{yo-y\)+ho{y'o+y{) 


(82) 


TESTING  IMPLEMENTATION  PARAMETERS  FOR  SECOND  ORDER 


A  test  equation  first  proposed  by  Prothero  and  Robinson  (Reference  27)  was 
extensively  utilized: 


y'  =  X(y-p(t))  +  p'(t),  y(t0)  =  y0.  (83) 

The  exact  solution  is 

y(t)  =  (>o  -  p{f o))exp(^(?  “  ?o))  +  P(0  •  (84) 


It  is  interesting  to  note  that  for  y0  =  p(t0)  the  solution  is  simply  p(t).  The  sine  function 


was  used  for  p(t)  and  the  interval  of  interest  is  t  e  [0,20],  Typical  values  for  X  were  -2  and 
-20. 


A  quantity  r  is  determined  for  each  step,  where  r  is  defined  as 


err 

est 


(85) 


Here  err  represents  the  local  error  as  described  in  Equation  40,  that  is,  the  solution  of  the 
initial  value  problem  in  Equation  1  beginning  at  a  DIMSIM  solution  point  (tn.I,yn.1)  with 
exact  local  solution  at  tn  given  by  y(tn;tn.1,yn_1),  and  with  calculated  solution  at  tn  of  yn. 
Thus  we  test  how  well  the  estimate  approximates  the  error  using  Equation  40, 
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errn=yn-y{tn’tn-\>yn-\)- 


It  should  be  noted  that  it  is  not  appropriate  to  restart  the  solution  in  seeking  to  obtain  the 
expression  for  err.  A  DIMSIM  changes  character  after  the  first  step,  and  this  must  be 
preserved  in  the  testing  process.  In  particular  the  external  stage  vector  is  not  recalculated 
using  a  starting  procedure,  otherwise  the  separate  error  estimate  for  an  intial  step  would  be 
required.  Testing  will  only  then  truly  reveal  the  behavior  of  the  error  estimator  in  its  use 
within  a  solver. 

The  results  of  a  test  may  be  indicated  both  graphically  and  using  a  number  of  statistics. 
The  quantity  est  is  used  to  denote  the  result  of  the  use  of  the  error  estimate  provided  for  the 
method.  A  graph  showing  both  err  and  est  together  is  sometimes  very  revealing,  but  it  is 
also  vital  to  look  at  the  largest  error  and  the  solution  range.  And  a  table  showing  the 
percentages  of  error  estimates  yielding  r  less  than  1%,  5%,  10%,  25%,  50%,  and  100% 
has  proven  to  be  of  great  value.  Of  course  it  is  expected  that  test  results  will  vary  greatly 
with  the  tolerance  used  and  also  with  the  problem  solved  and  even  with  the  problem 
parameters.  Thus  any  testing  of  implementation  parameters  must  be  considered  preliminary 
to  the  much  more  extensive  testing  required  of  a  full  solver. 

Two  tests  were  provided  for  the  quality  of  error  estimates.  Test  2  was  carried  out  as 
part  of  the  testing  of  the  complete  solver  on  a  number  of  different  problems  and  is 
described  with  results  in  the  Developing  Explicit  DIMSIM  ODE  Solvers  section.  Test  1 
was  first  described  by  Butcher  and  Jackiewicz  (Reference  3).  A  scheme  was  devised  for 
testing  the  effects  of  rapid  step  changes  on  error  estimation  accuracy.  An  initial  step  size  h0 
is  chosen,  along  with  a  ratio  p.  The  cyclic  pattern  of  step  sizes  h0,  ph0,  p2h0,  ph0,  h0,  hj 
p,  h0/p2,  h(/  p,  h0,  ...,  is  used  until  the  end  of  the  interval  of  interest  is  reached.  This 
should  be  considered  to  be  a  very  stringent  test,  especially  for  higher  values  of  p ,  and 
serves  as  an  excellent  preliminary  check.  The  behavior  even  with  rapid  step-size  changing 
is  quite  good,  as  shown  in  Tables  1  and  2,  and  Figure  1. 


TABLE  1.  Error  Estimate  Test  for  Second  Order,  Part  1. 


p 

%  r  <  0.01 

%r  <  0.05 

%r  <  0.10 

%r  <  0.  25 

%r  <  0.50 

%r  <  1.0 

1.25 

0.66 

9.27 

21.19 

49.67 

95.36 

99.34 

1.50 

3.29 

11.84 

23.68 

49.34 

91.45 

99.34 

1.75 

2.63 

6.58 

13.16 

34.87 

78.29 

99.34 

2.00 

0.65 

3.92 

9.80 

25.49 

67.32 

99.35 
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TABLE  2.  Error  Estimate  Test  for  Second  Order,  Part  2. 


p 

rmax 

tf 

max  err 

max 

err/Atj 

1.25 

0.2811 

3.2740 

20.0402 

0.0011 

0.0031 

6.6206 

20.0963 

0.0010 

1.75 

0.1719 

3.8194 

20.0335  | 

0.0010 

0.0402 

22.0300 

20.1128 

0.0010 

FIGURE  1 .  Graph  of  Error  Estimate  Test  for  Second  Order.  Estimate  is 
dotted  line,  local  error  is  solid  line. 


A  related  test  was  devised  for  interpolants.  The  interpolant  was  used  to  calculate  the 
solution  at  10  points  evenly  spaced  across  a  step.  The  error  at  the  end  of  the  step  was  used 
to  scale  the  error  at  the  interpolated  points.  Thus  a  ratio  no  greater  than  1  in  absolute  value 
would  indicate  that  the  error  at  the  interpolated  point  is  no  greater  than  the  error  at  the  end 
point.  Any  of  the  error  estimator  tests  described  above  might  be  used  to  generate  the 
integration  points,  but  the  rapid  step-size  changing  test  seemed  most  directly  useful. 

Figure  2  shows  the  behavior  of  the  Butcher- Jackiewicz-type  interpolant  (Equation  81) 
error  relative  to  the  behavior  of  the  solver  at  the  grid  points.  It  performs  very  well  with 
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only  a  few  of  the  thousands  of  interpolated  points  having  an  error  greater  than  the  error  at 
the  grid  points.  In  fact,  much  of  the  time  the  error  at  the  inteipolated  points  in  the  middle  of 
the  interval  is  less  than  the  end  point  error.  Certain  grid  points  seem  to  lie  close  to  a  zero  of 
the  error  function,  having  associated  errors  as  much  as  2  to  4  times  less  than  those  at  other 
points,  and  the  interpolated  values  in  the  middle  of  those  intervals  have  significantly  less 
accuracy  than  the  grid  points. 


RHO  =  1.25 


RHO  =  1.5 


-5 


-5 


FIGURE  2.  Test  Results  for  Order  2  Butcher- Jackie wicz-type  Interpolant. 


It  is  interesting  to  look  at  the  percentage  of  relative  errors  falling  in  various  ranges.  For 
-1500  steps  the  factors  k  by  which  the  interpolation  errors  were  greater  than  the  grid  point 
error  were  distributed  as  shown  in  Table  3. 
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TABLE  3.  Test  Results  for  Order  2  Butcher-Jackiewicz-Type  Interpolant. 


p 

k  >  10 

10>k>5 

5  >  k  >  2 

2>k>  1 

k<=  1 

1.25 

0.0% 

0.0% 

0.0% 

1.8% 

98.2% 

0.6% 

0.4% 

0.1% 

0.9% 

1.75 

0.0% 

0.1% 

0.8% 

0.4% 

97.7% 

0.3% 

0.2% 

0.4% 

1.5% 

97.6% 

The  Nordsieck  interpolant  was  also  tested  using  the  same  test  problem.  Figure  3 
shows  that  this  interpolant  has  errors  that  are  fairly  comparable  to  the  grid  point  errors. 
The  simple  second  order  Nordsieck  vector  was  used  without  the  third  order  continuity 
correction  term. 


RHO-1.25  RHO  =  1.5 


FIGURE  3.  Order  2  Nordsieck  Interpolant  Test. 


For  -1500  steps  the  factors  k  by  which  the  interpolation  errors  were  greater  than  the  grid 
point  error  were  distributed  as  shown  in  Table  4. 
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TABLE  4.  Order  2  Nordsieck  Interpolant  Test. 


k>  10  10>k>5  5  >  k  >  2  2>k>l  k<=l 


0.0% 

0.9% 

0.1% 

0.4% 


0.1% 

0.1% 

0.5% 

0.4% 


10.2% 

78.3% 

13.1% 

74.7% 

15.5% 

72.6% 

17.1% 

70.6% 

11.3% 

11.2% 

11.2% 

11.6% 


It  is  evident  that  for  this  method,  the  interpolant  of  Butcher-Jackiewicz  type  (Equation  81) 
provides  a  better  representation  of  the  solution  than  the  simple  2nd  order  Nordsieck 
interpolant. 

The  continuous  version  was  then  tried  with  the  graphical  results  shown  in  Figure  4. 


RHO  =  1.25 


RHO  =  1.5 


FIGURE  4.  Order  2  Continuous  Nordsieck  Interpolant  Test. 
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The  statistical  breakdown  obtained  by  continuous  Nordsieck  interpolant  test  is  shown  in 
Table  5. 


TABLE  5.  Order  2  Continuous  Nordsieck  Interpolant  Test. 


p 

k>  10 

10  >  k  >  5 

5  >  k  >  2 

2  >  k  >  1 

k<  =  1 

1.25 

0.0% 

0.0% 

1.2% 

68.1% 

30.7% 

0.6% 

0.4% 

0.5% 

67.8% 

30.7% 

1.75 

0.0% 

0.2% 

1.5% 

67.1% 

31.2% 

0.3% 

6.2% 

1.2% 

66.7% 

31.6% 

This  interpolant  is  less  accurate  here  than  the  Butcher-Jackiewicz-type  interpolant  (Equation 
81),  but  provides  substantially  better  accuracy  than  the  2nd  order  Nordsieck  interpolant  and 
could  be  readily  used  if  a  better  alternative  did  not  exist  for  this  case. 

The  testing  of  the  starting  method  concerns  the  accuracy  of  the  starting  vector,  the 
accuracy  of  the  first  step  error  estimate,  and  the  appropriateness  of  the  choice  for  initial  step 
size.  The  Prothero-Robinson  test  problem  is  used  here  with  A,  =  -2,  with  a  starting  point  at 
t0  on  the  exact  solution  trajectory  through  (0,  1). 

The  following  test  reveals  the  accuracy  of  the  error  estimate  and  of  the  initial  step  size 
selection  algorithm.  The  starting  point  att0=  1  is  used  throughout  to  avoid  the  special  case 
of  a  0  third  derivative  at  t0.  It  is  evident  that  accuracy  improves  until  round-off  error 
becomes  significant.  The  fact  that  the  error  term  is  O(lr)  is  evident  as  a  tightening  of  the 
tolerance  by  a  factor  of  1000  produces  changes  of  a  factor  of  10  in  h. 


TABLE  6.  Order  2  Starting  Procedure  Test  1. 


Tol 

h  (Nord) 

Loc  err 

Err  est 

Loc  err 

Err  est 

10'3 

1.95x10'’ 

-3.21X10-4 

-5.04x1  O'4 

3.09xl0'2 

-1.87xl0'6 

-1.99xl0'6 

io-6 

1.95x1 0'2 

-4.81xl0'7 

-5-OlxlO'7 

3.09x1 0'3 

-1.97x1  O'9 

-1.99xl0"9 

10-9 

1.95x1  O'3 

-4.98x10-’° 

-5.00x10-’° 

3.09X10-4 

-2.02xl0'12 

-2.00xl0'12 

10-12 

1.95x1  O'4 

-5.00xl0’3 

-5.00xl0"13 

3.09x1  O'5 

-1.89xl0‘5 

-1.97xl0'15 

10-15 

1.95x1  O'5  ! 

-5.55xl0"16 

-5.1  lxlO16 

3.09x1  O'6 

-l.llxlO'16 

-3.75xl0"’8 

Another  revealing  test  uses  starting  points  at  0  through  1,  evenly  spaced,  and  with  a 
fixed  tolerance  of  10  .  We  note  in  Table  7  that  the  use  of  the  Nordsieck- vector-based  error 
estimate  for  the  first  step  provides  a  much  more  efficient  start  than  the  method  described 
above  in  the  previous  chapter  as  outlined  in  Hairer,  Norsett  and  Wanner  (Reference  18). 
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TABLE  7.  Order  2  Starting  Procedure  Test  2. 


tf) 

■msH 

Loc  err 

Err  est 

Loc  err 

Err  est 

0 

l.lOxlO'03 

-4.99x1  O'10 

-5.00x10-'° 

1.36x1a04 

-9.38x1  O'13 

-9.39xl0"13 

0.2 

1.24x1  O'03 

-4.99x1  O'10 

-5.00x10'° 

1.59xlO"04 

-1.07xl0"12 

-1.07xl0'2 

0.4 

1.39xl0"°3 

-4.99x10"'° 

-5.00x10'° 

2.11x1a04 

-1.77xl0"12 

-1.78xl0’2 

0.6 

1.55xl0-°3 

-4.98xl0'10 

-5.00x10-'° 

2.55xl0"°4 

-2.24x1  O'12 

-2.24x10" 12 

0.8 

1.73xl0-°3 

-4.98x10-'° 

-5.00x10'° 

3.24x1a04 

-3.29xl0"12 

-3.29xl0"'2 

1 

1.95xl0"3 

-4.98x1  (f10 

-5.00x10'° 

3.09xia4 

-2.02xl0-‘2 

-2.00xl0‘2 

The  final  test  concerns  the  calculation  of  the  starting  vector.  Here  we  are  concerned  that 
the  smallest  possible  safe  value  for  h0  be  utilized  to  produce  the  maximum  accuracy  since 
no  DIMSIM  integration  step  will  actually  be  taken.  We  estimate  that  the  derivatives  will  all 

be  approximately  at  least  the  same  size  as  g(t0,y0)=min(||y0||,|/(r0,y0)|),  but  taking  the 

maximum  if  the  minimum  is  smaller  by  a  factor  of  eps,  the  machine  epsilon.  We  note  that 


quantities  h^y^k\t0)  will  be  calculated  in  terms  of  linear  combinations  of  functional  values 

and  derivatives  multiplied  by  h0,  with  the  most  difficult  term  to  calculate  corresponding  to 
k=p+l,  used  for  the  initial  step-size  selection.  We  would  like  to  have  at  least  three 
significant  digits  for  this  derivative,  and  we  expect  to  have  more  correct  digits  for  the  other 

terms.  Thus  we  would  like  to  have  g(t0,y0)h£+l  =  \05 eps  (103  instead  of  105  worked  well 


for  order  2  but  not  for  order  5,  this  is  a  bit  more  conservative  but  more  generally  useful). 


This  yields  an  optimal  choice  of  /iq  =  10p+1 


eps 


\p+i 


In  this  case  of  course,  p  =  2. 


.s(Wo)y 

We  test  this  by  calculating  the  relative  errors  in  the  second  and  third  derivatives  at  5 
uniformly  spaced  starting  points  between  0  and  1 ,  and  also  looking  at  the  maximum  norm 
of  the  error  in  the  starting  external  stage  vector.  Note  that  for  this  test  problem  the 
increment  was  quite  appropriate  (see  Table  8). 


TABLE  8.  Order  2  Starting  Procedure  Test  3. 


to 

h0 

y'Vo) 

Rel  err 

r(< o) 

Rel  err 

IM1 

0 

4.00 

jByaniBi 

-9.00 

1.69x10“* 

0.2 

■HEHlIiM 

2.48 

-6.34 

MBIIMi 

4.11  xlOJ4 

0.4 

■ililrliliia 

1.41 

-4.51 

■rmiia 

1.01  xIO17" 

0.6 

6.40x1  O'* 

3.64x10’' 

-3.23 

■irTnrM 

2.51  xl0"“* 

0.8 

■atvninM 

-2.31 

1 

WeWlliM 

twnam— 

-1.62 

Ffmv 
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IMPLEMENTATION  PARAMETERS  FOR  A  FIFTH  ORDER  EXPLICIT 
DIMSIM 

A  fifth  order  type  1  DIMSIM  derived  by  Butcher  and  Jackiewicz  (Reference  8)  was 
adapted  for  implementation  as  an  ODE  solver.  The  method  was  derived  using  numerical 
techniques  for  solving  the  large  system  of  nonlinear  equation  and,  to  ensure  favorable 
stability,  adding  a  requirement  that  the  form  of  the  stability  polynomial  reduce  to  the 
standard  5th  order  Runge-Kutta  stability  polynomial.  Here  again,  p  =  q  =  r  =  s  =  5,  U  =  I. 
The  method  coefficients  are  as  follows: 


c 


Ill 

4  2  4 


A  = 


0 

1.1765281703106688 

1.9805793463233191 

3.0532835108392395 

1.4193325269467698 


0 

0 

.4017181027378085 

.7349961462028882 

2.6534897473125331 


0 

0 

0 

.2672357626475791 

-2.2778532945468265 


0 

0 

0 

0 

1.1978905088172778 


0 

0 

0 

0 

0 


V  is  given  by  a  matrix  with  5  identical  rows,  where  each  row  is  vT,  and 


v  = 


'  -.2406956155386215 
1.2604945758471451 
-2.4812693924523267 
1.9199070083032958 
0.5415634238405073 


B  = 


'3.163023914364555736 

3.250176692142333514 

3.508528033848969459 

4.176223954923772036 

3.008220785304765914 


1.974339246050540681 

1.53197813493942957 

0.327374204184027625 

-2.61827171939311992 

2.11453341958770507 


-0.81042512005562813 

0.097908213277705203 

2.239060519232953536 

7.03900409877443037 

0.3572048837825639 


0.540922018802018401 

-0.422272425642426043 

-2.097452509375452155 

-5.2884360132659356 

-1.90425330857598604 


0.05507865419631683844' 

-0.4613800716699075171 

-0.936868983593822539 

-1.691097027371050162 

-0.64562655527099952 


and  was  calculated  using  Theorem  2  and  Mathematics  The  computation  of  W  is 
straightforward  using  the  first  relation  of  Equation  11, 
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1 


W 


z 

_2 


=  (I-zA)ecz  +  0(z6). 


This  yields,  equating  terms  of  the  same  degree, 


W  = 


e 


c-Ae 


or  evaluating  this  numerically, 


'1  o 

1  -0.926528170310669 
1  -1.882297449061128 
1  -3.305515419689707 
1  -1.992859488529754 


0 

0.03125 

0.02457047431554788 

-0.0361169178745116 

0.07713632883232162 


0  0  0 
0.002604166666666666  0.0001627604166666666  8.1380208333333310^-6 
0.00827964262277682  0.001558025774120291  0.0001950328608825182 

0.01393940010021236  0.005702129564105415  0.001161984318267553 

0.03157006827664394  -0.002014862315115681  -0.001959141967572072 


This  approach  was  utilized  in  calculating  the  rescaling  matrices  for  this  5th  order  DIMSIM. 

The  customary  form  of  V  is  utilized  with  a  first  row  of  vT  and  the  remaining  elements  0, 
while 


"3.163023914364555 

1.974339246050541 

-0.810425120055628 

0.5409220188020184 

0.055078654196316831 

0 

0 

0 

0 

1 

1 

_I6_ 

3 

12 

-16 

21 

3 

B  = 

44 

224 

152 

416 

140 

3 

3 

3 

3 

96 

-448 

768 

-576 

160 

256 

-1024 

1536 

-1024 

256 

was  computed  using  Mathematica. 

Interpolation  may  always  be  carried  out  using  the  Nordsieck  vector  or  the  continuous 
modification  of  the  Nordsieck  interpolant  as  described  above  in  the  Implementing 
DIMSIMs  section.  However,  the  Butcher-Jackiewicz  interpolant  has  some  desirable 
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features  when  it  exists,  and  it  was  also  calculated  for  this  method.  It  is  a  continuous 
interpolant  of  both  order  and  degree  5,  while  the  simple  Nordsieck  interpolant  is  also  of 
order  and  degree  5  but  not  continuous  at  the  grid  points,  and  the  continuous  Nordsieck 
inteipolant  is  of  degree  6  but  adds  the  additional  oscillation  typical  of  a  higher  degree 
polynomial  without  increasing  the  order  of  accuracy,  since  the  method  itself  is  only  of 
order  5.  The  interpolant  sought  is  of  the  form  Equation  34,  which  in  this  case  is 

>?(<„- 1  +e>h,)= Wo(0)f(i'w) + 

where  (30  and  y0  are  both  column  vectors  of  polynomials  of  degree  5.  We  write  them  in  the 


AoiO  +  Aoil^  +  A)1202  +  Poi3®3  +  A)14^4  +  PotfQ5 
A)20  +  Po2l®  +  Po22®2  +  Po23®3  +  P024Q4  +  P025®5 
Po  =  A)30  +  Po3l®  +  Po32®2  +  Po33®3  +  A)34#4  +  Po35Q5 
Po40  +  P 041®  +  Po42®2  +  Po43®3  +  Po44®4  +  Po45®5 
_A)50  +  Po51®  +  Po52®2  +  Po53®3  +  Po54®4  +  Po55®5  . 

and 

7oio  +  7oi  l®  +  7oi2^2  +  7oi3^3  +  7O1404  +  7oi5^5 
7020  +  7 021®  +  7 Oil®2  +  7o23®3  +  7 024® 4  +  7025®5 
7o  =  7030  +  7031^  +  7o32®2  +  7033®3  +  7034®4  +  7o35^5  * 

7040  +  7041^  +  7042®2  +  7o43^3  +  7 044® 4  +  7 045® 5 
_7050  +  7051^  +  7o52®2  +  7o53®3  +  7o54®4  +  7o55®5_ 

We  use  the  criterion  Equations  35  and  36  to  evaluate  the  unknown  coefficients.  That 
p0j0=  0  for  all  j  follows  immediately  from  the  condition  |30(0)  =  0.  From  Corollary  1  .ii  to 
Theorem  5  we  find  that  we  must  have  y0(0)  =  [l  0  0  0  0].  This  reduces  the 
continuity  conditions  of  Equation  36  and  to  the  following  relations: 

I  Ac ■)ij=h,b  i  =  1,2,..., 5 
7=1 

and 
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X Yoij  —  >  *  1)2,.. .,5. 
j= 0 


There  are  a  total  of  50  unknowns  remaining  and  10  linear  equations  have  already  been 
specified.  The  order  condition  in  Equation  35  provides  additional  linear  equations.  We 
have,  again, 


zp0(d)ecz  +  y0(e)w(z)  =  +  o(z6),  6  e  (0,1] . 


We  find  the  following  equations  for  terms  of  degree  0  through  5  in  z,  note  that  the  alphas 
are  the  elements  of  W  and  that  the  first  row  of  W  is  0  except  for  a  first  element  1 . 


DegreeO:  Y0(d)a0=y0(e)e  =  l. 


Degree  1:  p0(9)e  +  f,Yoi(e)a> :i  =  e- 
i=2 

5  5  q2 

Degre_e_2:  'LP0i(6)ci  +  lYoii^n  =  ~7T- 

i=2  i=2  2 


Degree  3:  X  Poi(d)^T  +  2  yo/C^s  =  7"- 

i=2  2  i'=2  O 


2  5 


9i 


Degree4:  X  Poi(9)^  +  X7o;(0)«*4  = 

i=2  6  i=2 


5  cf  5  _  0 


5 


Degree_5:  X  Poi(e)^7+  lYoii^iS  =  10n 
j=2  24  ;'=2  120 


Setting  terms  of  the  same  degree  in  q  in  each  of  these  equations  to  0  provides  5  equations 
for  each  degree,  for  a  total  of  30  additional  equations.  This  means  that  there  are  a  total  of 
40  equations  to  solve  with  50  unknowns,  leaving  10  free  parameters. 

The  equations  were  solved  using  Mathematica  using  a  process  of  symbolic  elimination 
of  variables.  The  decimal  representation  of  the  method  parameters  resulted  in  the 
appearance  of  some  inexact  floating  point  numbers.  When  the  coefficients  or  constant 
terms  were  l.Oe  -  16  or  less  they  were  understood  to  be  0  for  the  purpose  of  satisfying  the 
interpolant  defining  conditions.  It  might  be  desirable  to  also  try  solving  these  equations 
with  a  purely  numerical  solver.  The  following  interpolant,  one  of  a  family  of  equivalent 
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interpolants,  was  calculated  in  this  way,  with  free  parameters  identified  during  the  solution 
process  set  to  0  at  the  end.  It  should  be  noted  that  an  alternative  procedure  for  setting  free 
parameters  would  be  to  use  them  to  minimize  the  interpolant  error,  yielded  by  the  sixth 
degree  term  in  z  in  the  Taylor  expansion  of  Equation  35. 


Yo(e )  = 


1  +  6.6257606574898426 -2.49544672768745802 -1O.5O65833728538603  +  1O.4441633211653404  -  5.3085894936524990s]7" 

1 .260494575847 1450 
-2.48 1 2693924523270 

-2.9781439O48957180  +  4.898O5O913199O1302 

-2.4268419359889530  -  2.4026041 855 1 155602  +  10.506583372853860s  -10.444163321 1653404  +  5.3085894936524990s 


A>(o)= 


-14.438897214071 690  +  7. 15197332O6O746402  +  21. 8729345991698403  -  2O.678O835712726804  +  9.255096779931 650s  T 
-7.2412112973764320  +  4.5177685O937919802  +  1.O32326O5991517803  +  3.6654559741325980s 
6.0067479405540020  -  6.817173O6O6O96304 

-1.6721559444299640  -  3.594382030131 15303  +  9.869182843675O404  -  4.06172285031 19060s 
0.1622576215384850- O.26721155826OO82302 +1.6272657376669103  -  3.1876761861315804  +  1.7204430393825870s 


Although  there  seemed  to  be  only  10  free  parameters  based  on  the  number  of  equations  to 
solve,  some  of  the  equations  proved  to  be  dependent  as  they  became  effectively  0  before 
they  could  be  utilized  in  the  solution  process. 

The  error  coefficient  is  a  point  of  comparison  among  competing  methods  as  the 
coefficient  of  the  leading  discretization  error  term,  for  these  DIMSIMs  always  multiplying  a 

term  of  the  form  y^p+1\tn)h%+1 .  We  note,  first  of  all,  that  the  error  coefficient  of  any 
DIMSIM  is  given  by 


=  v 


v*=i  k\ 


Be p_ 
P]- 


\ 

J 


from  Theorem  6.  For  this  method  we  find  that 


<pP 


"0.000055 1 86676403623959 
0.000135380090032096182 
0.00045011570559905138 
0.00141313551999348304 
-0.0006734285188458039 
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and  so  the  error  coefficient  of  the  method  becomes  vT(pP  =  1.38889e-3.  However, 
interestingly  enough,  for  the  first  step  the  comparable  number  obtained  using  Theorem  8  is 
5.5 1 8667640362375e-05. 

Error  estimation  parameters  were  also  calculated.  The  initial  error  is  estimated  using 
vectors  P  and  y  where 


P  = 


'  0.1812875545254936 
0.310525433525683 6 
-0.3258295105315199 
0.1138267699195702 
-0.01051030018603054 


and 


7  = 


'-0.1351324309632461 

0 

0 

0 

0.1351324309632461 


These  are  calculated  using  the  conditions  of  Theorem  8,  which  yield  7  linear  equations  for 
the  10  parameters.  These  were  solved  using  Mathematica,  and  y2,  y3,  and  y4  were  free 
parameters  that  were  arbitrarily  set  to  0  to  produce  a  simple  form  minimizing  the  number  of 
operations  required  per  step. 

Subsequent  steps  utilize  different  vectors  p  and  y  where  these  are  determined  using 
Theorem  9.  Upon  examination,  the  first  condition  of  Equation  57  reduces  in  the  case  to  the 
second  and  so  there  are  again  7  linear  equations  in  10  unknowns,  solved  using 
Mathematica.  For  this  method  P  and  y  were  derived  to  be 


p  =  k\ 


-126.7321008977760 

-9.065592286799085 

9.939847085378934 

-1.234330257611486 

0 
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and 


y  =  k 


'  67.38322836913306  ' 

0 

-65.05966877019942 

0 

-2.323559598933673 


where 


_ 5^ _ 

1  + 5.7377413289581355 +  7.61378531497757652  +  2.929 17247339678653  +  . 053 12848737725 84 154  ' 


Observe  that  k  has  four  poles,  at  8  =  -52.4401,  8  =  -1.44533,  8  =  -1,  and  8=  -0.248337, 
and  these  negative  poles  are  not  a  problem  since  8  must  be  positive  in  any  ODE  solver. 

However,  it  should  be  noted  that  there  was  a  choice  here  of  free  parameters,  which 
were  all  arbitrarily  set  to  0  at  the  end  of  the  calculation.  These  were  p3,  y2  and  y4.  Other 
choices  of  free  parameters  were  tried  first  and  actually  yielded  poles  for  small  positive  8. 
In  particular,  a  choice  of  P5,  y3  and  y4  yielded  a  pole  at  8  =  0.93028;  a  choice  of  Ps,  y2 
and  y3  yielded  a  pole  at  8  =  0.0752006;  and  a  choice  of  y2,  y3  and  y4  yielded  a  pole  at 
8  =  0.568585.  Use  of  free  parameters  to  enforce  the  condition  pTe=0  also  resulted  in  a 
small  positive  pole  and  so  this  condition  was  not  enforced.  Testing  showed  that  so  long  as 
the  pole  was  avoided,  similar  results  were  obtained  for  all  choices  of  free  parameters. 
Clearly  it  is  undesirable  for  a  pole  to  appear  in  the  error  estimation  expression  in  the 
neighborhood  of  likely  choices  for  successive  step  size  ratios,  and  it  is  evident  that  these 
poles  are  unfortunate  artifacts  that  do  not  reflect  die  actual  errors  produced  through  step- 
size  changes.  This  problem  was  earlier  observed  by  Jackiewicz  and  Zennaro  (Reference 
28)  and  by  Bellen,  Jackiewicz  and  Zennaro  (Reference  29)  in  the  context  of  two-step 
Runge-Kutta  methods.  The  success  here  in  using  the  parameter  freedom  to  move  these 
poles  to  the  negative  real  axis  is  an  encouraging  result  that  may  have  direct  relevance  to 
other  ODE  methods. 

For  our  explicit  fifth  order  solver,  we  need  a  total  of  5  equations  to  obtain  for  our  5 
unknowns  the  derivatives  2  through  6.  We  use  an  order  6  explicit  Runge-Kutta  solver  to 
obtain  an  approximation  y,  to  y(t0+h0)  accurate  to  O(h07).  Then  we  have: 


>1  =  >o + Vo +^ty'o +^y'o  +  24y^ + + m^o6) + °{fi) 

y{  =  y'o+ Vo + t‘3;o"+  x  ^o4)  +Tbo5) + +  °(^  )• 


76 


NAWCWPNS  TP  8340 


Here  y{  =  f(t0  +  h^,y]).  Similarly  we  obtain  two  more  equations  by  computing  an 
approximation  y2  to  y(t0+h0/2)  accurate  to  O(h07).  This  yields 


yi=yo+bo+iyo*iyS'+myo)+moyo)+mioyo)+o(ii) 


v'  —  v'  +  v"'+  —  v(4)  +  v<5>  +  —  + 

y2  —  Jo  2  -^0  *  °  ***  '  4<?  yes  ^  ao/i  rn  *  ao/in  yn  ^ 


±  v"'  .  &  Vv  ■ 

8  ^0  +  48  ^0 


384  %  +  3840  ?0 


o(4\ 


Note  that  the  6th  derivative  term  must  be  included  to  obtain  sufficient  accuracy,  even  if  an  a 
priori  error  estimate  is  not  sought.  This  means  that  one  more  equation  must  be  used,  and 
this  may  come  from  computing  an  approximation  y4  to  y(t0+h0/4)  accurate  to  O(h07).  This 
final  equation  is 


hn 


1,2 


Eliminating  variables  we  solve  these  five  simultaneous  equations  and  obtain 


y'o  =  A- (-567 y0  -  25yi  -  432y2  +  1024y4  +  ^(-90^  +  3y{  +  72y'2))  +  o(^), 

9/zq 

y'o  =  4(l836y0  + 148^  +  21 12 y2  -  4096y4  +  h0(222y'0  - 1 8y{  -  384^))  +  C>(/^) 

k) 

y{0]  =  ^r(-1323yo  -  1836y2  +  3328y4  +  ^(-lUy'o  +  2ly{  +  378^))  +  o(h%) 

3/zo 

y(05)  =  ^(378 y0  +  70^  +  576y2  - 1024 y4  +  hQ(39y'0  -  9y{  - 132^))  +  o(h%) 

"0 

yi6)  =  ^(-90y0  -  22y,  -144y2+  2S6yt  +  ko(-9yi  +  3y{  +  36^))  +  O^). 


The  error  in  the  first  step  is  given  by  Theorem  8,  with  h,  the  first  DIMSIM  step  size,  to  be 


77 


NAWCWPNS  TP  8340 


Ite i  = 


'l  5  cp 

--  X  jBi/  — 

6!  A  7  5!y 


^(fo)+c>(^i7) 


=  (7.063894579663840e  -  2)<56(-90y0  -  22^  -  144y2  +  256y4  +  /zo(-9y6  +  3 y{  +  36y2)) 


If  we  set  this  error  to  half  the  tolerance  %  and  use  a  norm  to  include  the  possibility  of 
systems  of  equations,  we  then  may  write  a  conservative  but  hopefully  accurate  choice  for 
initial  step  size  to  be  8h0,  where 


^2(7.063894579663840e  -  02)||-90y0  -  22^  -  144y2  +  256y4  +  ^(-9^  +  3y(  +  36y£)|  J 


TESTING  IMPLEMENTATION  PARAMETERS  FOR  FIFTH  ORDER 

a.  Variable  Step  Error  Estimate 

The  error  estimate  was  tested  in  the  same  way  for  fifth  order  as  for  second  order.  For 
the  same  time  steps  the  following  results  were  obtained  in  Figure  5. 
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Error  Estimate 
_  Local  Error 


FIGURE  5.  Order  5  Error  Estimation  Test. 


The  statistics  obtained  during  the  fifth  order  error  stimation  test  are  contained  in  Tables  9 
and  10. 


TABLE  9.  Order  5  Error  Estimation  Test  Part  1. 


p 

%  rc.Ol 

%r<.05 

%r<.  10 

%r<.25 

%r<.50 

%r<1.0 

1.25 

0.66 

3.31 

6.62 

17.88 

52.32 

88.08 

1.50 

1.32 

8.55 

15.13 

29.61 

54.61 

1.75 

0.66 

3.95 

8.55 

16.45 

38.16 

90.79 

0.00 

0.65 

5.88 

11.11 

25.49 

77.77 

79 


NAWCWPNS  TP  8340 


TABLE  10.  Order  5  Error  Estimation  Test  Part  2. 


p 

rmin 

rmax 

_ if _ 

Max  err 

Max  err/At5 

1.25 

2.74e-2 

16.85 

20.0402 

1.54e-8 

3.78e-4 

1.50 

2.99e-2 

96.67 

20.0963 

1.59e-8 

4.02e-4 

1.75 

1.76e-2 

12.96 

20.0335 

1.50e-8 

3.67e-4 

2.00 

5.84e-2 

26.48 

20.1128 

3.21e-5 

7.80e-l 

The  error  here  is  significantly  less  than  what  appeared  with  order  2,  but  the  problem  for 
p  =  2  with  larger  errors  and  change  in  solution  order  is  troubling.  A  separate  test  was  run 
with  smaller  step  size,  a  difference  of  a  factor  of  two.  The  results  can  be  seen  in  Figure  6, 
and  Tables  11  and  12. 


x10-i°  RHO  =  1.25 


x10-10  RHO  =  1.75 


XlO 


r10  RHO  =  1.5 


Error  Estimate 
_  Local  Error 


FIGURE  6.  Order  5  Error  Estimation  Test,  Smaller  Step  Size. 
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TABLE  11.  Order  5  Error  Test,  Smaller  Step  Size,  Part  1. 


p 

%  rc.01 

%r<.05 

%r<.10 

%r<.25 

%r<.50 

%r<1.0 

1.25 

0.66 

4.98 

9.97 

27.24 

61.46 

80.73 

1.65 

9.27 

15.56 

29.14 

58.28 

1.75 

0.33 

4.30 

7.28 

18.87 

40.40 

92.72 

2.00 

0.66 

2.97 

5.61 

13.20 

33.66 

74.59 

TABLE  12.  Order  5  Error  Test,  Smaller  Step  Size,  Part  2. 


P 

rmin 

rmax 

U 

Max  err 

Max  err/At5 

1.25 

6.55e-4 

490.6 

20.0200 

4.96e-10 

3.83e-4 

1.50 

2.83e-2 

73.26 

20.0483 

5.10e-10 

3.99e-4 

1.75 

1.31e-2 

41.03 

20.0165 

4.88e-10 

3.76e-4 

2.00 

7.15e-3 

35.51 

20.0556 

5.11e-10 

3.92e-4 

We  can  see  that  in  this  case  the  results  are  consistent  for  all  p. 
b.  Interpolants 

Interpolation  testing  was  similar  to  the  procedures  used  for  the  second  order 
method,  but  with  more  points  and  over  a  shorter  interval.  The  results  for  the  Butcher- 
Jackiewicz  interpolant  are  seen  in  Figure  7. 
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RHO  =  1 .25  RHO  =  1.5 


FIGURE  7.  Test  Results  for  Order  5  Butcher-Jackiewicz-Type  Interpolant. 


The  statistics  obtained  during  the  fifth  order  interpolant  test  are  contained  in  Table  13. 


TABLE  13.  Order  5  Butcher-Jackiewicz  Interpolant  Test. 


p 

k>  10 

10  >  k  >  5 

5  >  k  >  2 

2  >k>  1 

k<=  1 

1.25 

0.0% 

0.0% 

0.4% 

8.8% 

1.50 

0.2% 

0.2% 

3.5% 

20.9% 

75.1% 

1.75 

2.4% 

2.0% 

7.5% 

28.1% 

59.9% 

2.00 

1.2% 

0.8% 

6.5% 

32.6% 

58.8% 

Evidently  some  significant  degradation  occurs  at  this  higher  order  with  increasingly  rapid 
step  changing.  A  similar  phenomenon  was  also  observed  in  a  study  of  related  general  linear 
methods  by  Enenkel  and  Jackson  (Reference  30)  and  Enenkel  (Reference  21).  It  should  be 
again  noted  that  no  error  minimization  was  carried  out  using  the  free  parameters. 
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The  simple  fifth  order  Nordsieck  interpolant  yielded  the  following  in  Figure  8. 


RHO  =  1.25  RHO  =  1.5 


5 

0 

_5 

- 1  - ““ ll 

5  1 

RHO  =  1.75 

c 

5 

5  10 

RHO  =  2 

1 1  II  ■ 

5  1 

0 

-5 

0  c 

inf 

5  10 

FIGURE  8.  Order  5  Nordsieck  Interpolant. 


The  statistics  obtained  during  the  fifth  order  Nordsieck  Interpolant  test  are  contained  in 
Table  14. 


TABLE  14.  Order  5  Nordsieck  Interpolant. 


P 

k>  10 

10  >  k  >  5 

5  >  k  >  2 

2  >k>  1 

k<=  1 

1.25 

0.0% 

0.0% 

0.0% 

88.5% 

11.5% 

0.0% 

0.0% 

0.0% 

86.5% 

13.5% 

1.75 

0.0% 

0.0% 

0.1% 

84.5% 

15.4% 

0.0% 

0.0% 

0.0% 

85.5% 

14.5% 

This  is  clearly  highly  accurate  for  this  range  of  step  sizes. 
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The  continuous  Nordsieck  interpolant  displayed  even  greater  accuracy  in  Figure  9. 


RHO  =  1.75 


RHO  =  2 


FIGURE  9.  Order  5  Continuous  Nordsieck  Interpolant. 


The  statistics  obtained  during  the  fifth  order  continuous  Nordsieck  interpolant  test  are 
included  in  Table  15. 


TABLE  15.  Order  5  Continuous  Nordsieck  Interpolant. 


p 

k>  10 

10  >  k  >  5 

5  >  k  >  2 

2  >  k  >  1 

k<=  1 

1.25 

0.0% 

0.0% 

0.0% 

1.1%  ”1 

98.9% 

1.50 

0.0% 

0.0% 

0.0% 

1.3% 

98.7% 

1.75 

0.0% 

0.0% 

0.0% 

1.4% 

98.6% 

2.00 

0.0% 

0.0% 

0.0% 

1.3% 

98.7% 
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This  must  be  considered  highly  desirable,  perhaps  with  more  accuracy  than  is  typically 
needed. 

It  then  became  of  interest  to  test  the  results  for  the  larger  step  sizes  used  for  2nd  order. 
The  results  for  the  Butcher-Jackiewicz  interpolant  were  as  follows  in  Figure  10. 


RHO  =  1 .25  RHO  =  1.5 


FIGURE  10.  Order  5  Butcher-Jackiewicz  Interpolant  Test,  Longer  Step  Size. 

The  statistics  obtained  during  the  fifth  order  Butcher-Jackiewicz  inteipolant  test  are 
contained  in  Table  16. 


TABLE  16.  Order  5  Butcher-Jackiewicz  Interpolant  Test,  Longer  Step  Size. 


p 

k>  10 

10  >  k  >  5 

5  >  k  >  2 

2  >  k  >  1 

k  <  =  1 

1.25 

3.9% 

3.4% 

13.4% 

26.4% 

52.9% 

3.4% 

3.8% 

13.9% 

33.9% 

45.1% 

1.75 

3.6% 

3.4% 

6.8% 

24.7% 

61.5% 

2.5% 

2.5% 

4.5% 

21.9% 

68.6% 
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For  the  fifth  order  Nordsieck  interpolant  the  results  are  shown  in  Figure  11. 

RHO  =  1.25  RHO  =  1.5 

r _ : _  r _ 


0  10  20  0  10  20 
RHO  =  1.75  RHO  =  2 


0  10  20  0  10  20 


FIGURE  11.  Nordsieck  Interpolant  Test,  Longer  Step  Size. 


The  statistics  obtained  during  the  Nordsieck  interpolant  test  are  contained  in  Table  17. 


TABLE  17.  Nordsieck  Interpolant  Test,  Longer  Step  Size... 

10  >  k  >  5  |  5  >  k  >  2 

0.0%  I  (B% 

0.0%  0.3% 

0.0%  0.3% 

0.0%  0.0% 


For  the  continuous  Nordsieck  interpolant  we  have  Figure  12  and  Table  18. 


k>  10 


0.0% 


0.0% 


0.0% 


0.0% 


2  >  k  >  1 


80.7% 


83.0% 


86.8% 


88.8% 
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RHO  =  1.25 


p 

■ 

Iff 

li 

0 

RHO 

10  20 

=  1.75 

11 

l|L  ILy 

RHO  =  1.5 


npni 

inn 

0  1 

RHO 

0  20 

=  2 

[■ 

pm 

0  20 

FIGURE  12.  Order  5  Continuous  Nordsieck  Interpolant  Test,  Longer  Step  Size. 


TABLE  18.  Order  5  Continuous  Nordsieck  Interpolant  Test,  Longer  Step  Size. 


P 

f>  10 

10  >  f  >  5 

5  >  f  >  2 

2  >  f  >  1 

f<=  1 

1.25 

0.5% 

0.1% 

0.0% 

25.8% 

73.6% 

0.5% 

0.1% 

0.3% 

26.7% 

72.4% 

1.75 

0.0% 

0.0% 

0.0% 

26.5% 

73.5% 

0.0% 

0.0% 

0.0% 

26.3% 

73.7% 

Both  Nordsieck  interpolants  seem  to  be  highly  acceptable,  with  the  simple  Nordsieck 
interpolant  producing  fewer  large  errors  with  k  >  5,  while  the  continuous  Nordsieck 
interpolant  otherwise  generating  higher  accuracy.  The  continuous  Nordsieck  interpolant 
seems  to  consistently  produce  k  <  1,  while  the  simple  Nordsieck  interpolant  produces 
many  values  of  k  slightly  greater  than  1 . 
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The  starting  procedure  testing  was  similar  to  that  used  for  the  second  order  method. 
The  selection  of  h0  proved  to  be  a  bit  delicate  and,  as  alluded  to  above,  the  formula 

ghfi  -  105  eps  proved  to  be  necessary  to  provide  a  suitable  value.  If  the  eps  factor  was 
reduced  by  a  factor  of  10  the  step  was  too  small  to  prevent  serious  round-off  errors  from 
occuring,  while  if  the  factor  was  increased  by  10  there  was  a  loss  of  local  truncation 
accuracy.  This  is  not  because  of  the  breakdown  in  the  assumption  that  the  derivatives  are 
of  the  same  size  as  the  initial  value  and  the  corresponding  derivative.  In  fact  in  this  test  6th 
derivatives  grew  to  be  as  large  as  64  while  first  derivatives  and  functional  values  were 
around  0.3  to  1,  but  this  would  tend  to  reduce  rather  than  increase  the  size  of  h0.  A  better 
answer  is  that  the  amount  of  cancellation  in  the  calculation  of  the  Nordsieck  vector 
components  requires  that  more  digits  exist  in  the  answer. 

The  testing  of  the  starting  method  concerns  the  accuracy  of  the  starting  vector,  the 
accuracy  of  the  first  step  error  estimate,  and  the  appropriateness  of  the  choice  for  initial  step 
size.  The  Prothero-Robinson  test  problem  is  again  used  here  with  X  =  -2,  with  a  starting 
point  at  ^  on  the  exact  solution  trajectory  through  (0,1). 

The  following  test  (Table  19)  reveals  the  accuracy  of  the  error  estimate  and  of  the  initial 
step  size  selection  algorithm.  The  starting  point  at  t0  =  1  is  used  here  throughout.  It  is 
evident  that  accuracy  improves  until  round-off  error  becomes  significant. 


TABLE  19.  Order  5  Starting  Procedure  Test  1. 


to 

■m 

Loc  err 

Err  est 

Loc  err 

Err  est 

10'3 

1.03 

-5.22x1  O'5 

6.81X10-4 

1.76x10'* 

8.05xl0'9 

1.34xl0'8 

10'6 

3.26x10* 

2.12xl0'7 

5.77xl0'7 

5.55xl0'2 

1. 09x10* 1 

1.26x10'** 

10'9 

1. 03x10'* 

3.96x10'*° 

5.32x10'*° 

1.76x1  O'2 

1.90xl0'14 

1.12xl0‘4 

10'*2 

3.26xl0'2 

5.09xl0'13 

5.00xl0'*3 

5.55xl0'3 

5.55xl0'16 

-2.39xl0'*7 

10-15 

1.03x1  O'2 

1.89xl0'15 

3.26xl0'16 

1.76x1 0'3 

5.55xl0'16 

-3.33xl018 

Another  revealing  test  uses  starting  points  at  0  through  1,  evenly  spaced,  and  with  a 
fixed  tolerance  of  10  .  At  this  tolerance  and  at  this  order  the  error  estimate  is  not  quite  as 
accurate.  We  note  again  that  the  use  of  the  Nordsieck-vector-based  error  estimate  for  the 
first  step  provides  a  much  more  efficient  start  than  the  method  outlined  in  Hairer,  Norsett 
and  Wanner  (Table  20). 


TABLE  20.  Order  5  Starting  Procedure  Test  2. 


to 

■Bl 

Loc  err 

Err  est 

Loc  err 

Err  est 

0 

7.22x1  O'2 

4.23x10'*° 

5.16x10'*° 

1.17xl0'2 

8.77xl0'15 

8.71xl0'15 

0.2 

7.75xl0'2 

4.29x10'*° 

5.26x10'*° 

1.26x1  O'2 

1.39xl0'14 

9.59xl0'5 

0.4 

8.32xl0'2 

4.16x10'*° 

5.28x10'*° 

1.45xl0'2 

1.20xl0'14 

1.48xl0'*4 

0.6 

8.91xl0'2 

4.12x10'*° 

5.31x10'*° 

1.60x1  O'2 

1.65xl0'14 

1.71xl0'14 

0.8 

9.54xl0'2 

4.14x10'*° 

5.28x10'*° 

1.80x1  O'2 

2.93xl0'*4 

2.36xl0'*4 

1 

1.03x10* 

3.97x10'*° 

5.32x10'*° 

1.76x1  O'2 

1.90xl0'14 

1.12xl0'12 
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The  final  test  concerns  the  calculation  of  the  starting  vector.  Here  we  are  concerned  that 
the  smallest  possible  safe  value  for  h0  be  utilized  to  produce  the  maximum  accuracy  since 
no  DIMSIM  integration  step  will  actually  be  taken  until  this  is  rescaled.  As  indicated 


above,  we  use  a  choice  of  Hq  -  lO^1 


eps 


V+i 


s(Wo). 


,  with  p  =  5.  We  test  this  by  calculating 


the  relative  errors  in  the  fifth  and  sixth  derivatives  (noting  that  higher  derivatives  are 
computed  less  accurately)  at  5  uniformly  spaced  starting  points  between  0  and  1.  The  error 
in  the  starting  vector  is  also  calculated  (Table  21). 


TABLE  21.  Order  5  Starting  Procedure  Test  3. 


to 

h0 

><5,('o) 

Rel  err 

Rel  err 

IBB 

0 

bhiiewM 

E— 

4.51xl0“3 

WSMMKSMM, 

ililUBM 

1.47x1 0'4 

mxuvkittm 

1.72  xlO'14 

0.4 

3.16xl0'z 

6.31xl0‘4 

tianiiM 

4.48  xlO'IJ 

0.6 

TTTTTTTMi 

-8.81 

BBE IBM 

■BdrflUiM 

beeehshi 

0.8 

WTHTf 

-5.76 

iti-nna 

WESSSSM 

1.48  xlO 

1 

KESsEIBMi 

-3.79 

EMM 

7.82 

■awaaiaa 

2.48  xlO"4 

THE  DIMEX  FAMILY  OF  EXPLICIT  DIMSIM  ODE  SOLVERS 


DESCRIPTION  AND  USE  OF  DIMEX  FAMILY  OF  SOLVERS 

The  second  order  type  1  DIMSIM  explicit  variable  step-size  ordinary  differential 
equation-initial  value  problem  solvers  DEMEXx  are  written  as  a  collection  of  double 
precision  FORTRAN  77  subroutines  called  by  a  driver  dimxl  that  provides  the  interface  to 
the  user’s  calling  program.  Here  x  denotes  the  method  order.  These  are  research  codes, 
designed  primarily  for  use  in  waveform  relaxation  in  a  form  that  should  make  them  suitable 
for  extension  to  higher  orders,  as  well  as  for  use  in  solving  standard  ODE  initial  value 
problems.  A  number  of  elements  remain  to  be  determined  in  developing  mature  production 
codes,  and  optimization  for  speed  and  memory  usage  have  not  been  performed.  The  codes 
use  the  Type  1  schemes  along  with  the  coefficients  described  in  the  previous  chapter. 
Advantage  is  taken  of  the  FASAL  property  referred  to  in  the  Implementing  DIMSIMs 
section  whenever  use  of  the  same  step  size  is  continued  for  another  step  for  both  order  2 
and  order  5,  and  for  all  steps  after  the  first  for  second  order.  In  this  section  we  discuss 
other  significant  elements  of  the  software  design  and  use.  The  complete  code  for  DIMEX5 
appears  in  the  appendix. 
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a.  Calling  the  Solvers 

The  user  must  issue  the  command  (x  is  the  order  number,  currently  2  or  5) 
call  dimxl(X0,Y0,X,Y,F,NEqn,XF,H, ATol,RTol, Starting, XOut, 

1  NCalls,NMisses,Diag,ErrEst,Work) 

Parameter  definitions: 

X0:  (real*8)  Initial  value  of  independent  variable  at  start  of  integration,  input 
Y0:  (real *8  array)  Initial  value  of  solution  at  start  of  integration,  dimensioned 
NEqn,  input 

X:  (real*8)  Desired  output  point  if  Diag  is  .false.,  otherwise  set  to  same  as  XF, 
input 

Y:  (real*8)  Computed  solution  at  X  if  Diag  is  .false.,  otherwise  at  XOut,  output 
only 

F:  (subroutine)  Provides  derivative  function  f  for  ODE  y’=f(x,y).  Subroutine  must 
be  defined  with  parameters  F(X,Y,YP,NEqn),  where  YP  is  the  derivative  at 
(X,Y)  and  NEqn  is  as  defined  below.  F  must  be  declared  external  in  the 
calling  program. 

NEqn:  (integer)  Number  of  equations  in  ODE  system,  input 
XF:  (real* 8)  Intended  termination  point  for  integration  process,  input 
H:  (real*  8)  Step  size  used.  Output  only,  not  to  be  defined  or  changed  by  the  user 
ATol:  (real*8)  Absolute  tolerance  for  local  truncation  error,  used  in  controlling 
integration  according  to  formula  described  below,  input 
RTol:  (real*8)  Relative  tolerance  for  local  truncation  error,  used  in  controlling 
integration  according  to  formula  ErrEst<ATol+RTol*  flTl^ ,  input 
Starting:  (logical)  Initially  set  to  .true.,  thereafter  not  to  be  changed  by  the  user 
XOut::  (real*8)  End  point  of  last  integration  step  performed,  output  only,  not  to  be 
changed  by  the  user 

NCalls:  (integer)  Number  of  function  calls  since  beginning  of  integration.  Output 
only,  not  to  be  changed  by  the  user. 

NMisses:  (integer)  Number  of  tolerance  misses  since  beginning  of  integration. 

Output  only,  not  to  be  changed  by  the  user. 

Diag:  (logical)  Set  to  true  if  only  one  integration  step  at  a  time  is  desired,  which  is 
the  preferred  mode  for  study  of  integrator  behavior,  and  false  if  it  is  desired 
to  integrate  until  interpolated  output  at  X  may  be  obtained 
ErrEst:  (real*8)  Maximum  norm  of  estimated  error  in  last  integration  step  taken, 
output  only 

Work:  (real*8)  Work  array  which  must  be  dimensioned  at  least  2+24*NEqn  for 
2nd  order  and  2+53*NEqn  for  5th  order. 
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b.  Calling  tree  for  DIMEX5 


dim51 

I 

driver 


start  solver  NewStep  Nordlnterp  rescale 

I  \  I 

F  RK6  F 

I 

F 

-A  block  data  subprogram  called  matrices  is  also  used. 

FIGURE  13.  Calling  Tree  for  DIMEX5. 


c.  Description  and  parameters  for  subprograms  of  DIMEX5 

subroutine  dim51:  See  Section  a  above  for  discussion  of  parameters.  This  subroutine 
simply  serves  the  purpose  of  hiding  from  the  user  the  details  of  the  numerous  arrays 
used  by  DIMEX5.  A  single  large  work  array  is  split  appropriately  as  required  by  the 
driver. 

subroutine  driver(XO,YO,X,Y,F,NEqn,XF,H,ATol,RTol,Starting,XOut, 

1  NCalls,NMisses,Diag?ErrEst,H01d,YIter,FStage,YIterP, 

2  X01d,YP,Yl,Y2,Y4,YlP,Y2P,YStage, 

3  YIterN,  YIterS  ,FS  tageN,  Y  6Der  ,FStageP,RKW  ork,N  ordV  ec) 

driver  provides  the  essential  logic  to  control  the  integration  process  and  calls  the 
subprograms  needed  to  carry  out  the  necessary  calculations.  When  Starting  is  set  to 
.true.,  start  is  called  to  obtain  an  initial  Nordsieck  vector,  step  size,  and  external 
stage  vector.  Solver  then  is  called  to  carry  out  the  first  integration  step  and  the  error 
estimate  is  checked  to  ensure  the  result  is  within  tolerance.  If  not,  the  step  size  is 
halved  and  Nordsieck  vector  is  rescaled  to  produce  a  new  starting  external  stage 
vector.  This  process  is  repeated  until  successful.  Starting  is  then  set  to  .false,  and 
the  integration  process  proceeds  as  NewStep  calculates  a  new  step  size  to  use.  If 
Diag  is  .true,  driver  returns  dim51,  which  returns  to  be  called  again  by  the  main 
program,  otherwise  at  each  iteration  either  XOut  is  at  least  X,  in  which  case 
Nordlnterp  produces  an  interpolated  solution  value  at  X  and  a  return  is  made  to  the 
main  program  for  another  call,  or  an  integration  step  is  carried  out.  An  integration 
step  begins  by  saving  copies  of  the  current  and  previous  external  stage  vectors. 
Then  rescale  is  called  to  rescale  the  current  external  stage  vector,  a  DIMSIM  step  is 
performed  using  a  call  to  solver,  and  ErrEst  is  checked  to  be  sure  the  error  is  within 
tolerance.  If  not,  a  tolerance  miss  is  counted  and  the  step  size  is  halved.  Rescaling 
again  takes  place  using  the  saved  values,  and  this  process  continues  until  the  step  is 
within  tolerance  or  MaxTries  is  exceeded,  generating  an  error  message  and  STOP. 
For  normal  execution  final  XOut  is  XF. 
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XO,YO,X,Y,F,NEqn,XF,H,ATol,RTol,Starting,XOut,NCalls,NMisses,Diag,ErrEst:  As 
in  dim51,  see  Section  a  above. 

HOld:  (real *8)  Step  size  used  for  previous  successful  integration  step 

YIter:  (real  *8  array)  External  stage  vector  for  current  step,  dimensioned  (NEqn,5) 

FStage:  (real *8  array)  Vector  of  derivative  values  computed  from  internal  stage  values  for 
most  recent  successful  step,  dimensioned  (NEqn,5). 

YIterP:  (real*  8  array)  External  stage  vector  for  previous  successful  step,  dimensioned 
(NEqn,5) 

XOld:  (real*8)  Previous  successful  step  point 
YP:  (real*8  array)  Derivative  at  X0,  =f(X0,Y0).  Dimensioned  NEqn. 
Y1,Y2,Y4,Y1P,Y2P,  Y6Der:  (real*8  arrays)  Storage  needed  for  subroutine  start.  Each 
dimensioned  NEqn. 

YIterN,YIterS,  FStageN:  (real*8  arrays)  Storage  needed  for  temporary  storage  of  current 
and  previous  external  stage  vectors  and  current  vector  of  derivatives  at  internal  stage 
points  until  it  is  verified  that  step  was  successful,  dimensioned  (NEqn,5). 

YStage:  (real*8  array)  Storage  needed  for  subroutine  solver,  dimensioned  NEqn. 

FStageP:  (real*8  array)  Previous  vector  of  derivative  values  for  last  internal  stage,  used 
again  for  first  internal  stage  of  new  step  if  step  size  does  not  change.  Dimensioned 
NEqn) 

RKWork:  (real*8  array)  Storage  needed  for  subroutine  start,  dimensioned  8*NEqn. 
NordVec:  (real*8  array)  At  first  step,  Nordsieck  vector  at  X0.  Subsequently,  Nordsieck 
vector  at  XOut.  Dimensioned  (NEqn, 6) 

subroutine  solver(X0,NEqn,YIter,H,HOld,F,ErrEst,FStage,YIterP,Starting,YStage, 

1  FStageP,KCalls) 

solver  applies  the  basic  DIMSIM  algorithm  to  compute  a  new  external  stage  vector 
and  an  error  estimate.  Different  coefficients  are  used  in  computing  die  error 
estimate  if  Starting  is  .true. 

X0:  (real*8)  Starting  point  for  this  integration  step 

NEqn, H,H01d,F,ErrEst, Starting:  See  description  of  driver,  above. 

FStage:  (real*8  array)  Upon  output,  set  to  vector  of  derivatives  at  internal  stage  points  for 
current  step,  dimensioned  (NEqn, 5) 

YIter:  (real*8)  Upon  input,  set  to  rescaled  previous  external  stage  vector.  Upon 
output,  set  to  current  external  stage  vector.  Dimensioned  (NEqn,5) 

YIterP:  (real*  8  array)  Output  only,  set  to  input  value  of  YIter.  Dimensioned  (NEqn, 5) 
YStage:  (real*8  array)  Temporary  storage  needed  for  each  internal  stage  vector  in  turn 
during  calculation.  Dimensioned  NEqn. 

FStageP:  (real*8  array)  Input,  see  driver  description  above.  Dimensioned  NEqn. 

KCalls:  (integer)  Number  of  calls  to  derivative  function  subroutine  for  this  solver  step, 
output 

subroutine  start(YIterP,F,X0,Y0,NEqn,H,YP,Tol,Yl,YlP,Y6Der,Y2,Y2P,Y4, 

1  YNord,RKW  ork) 

Start  first  calculates  the  derivative  at  the  initial  point  X0.  A  scaling  factor  g  to 
estimate  the  size  of  the  initial  Nordsieck  vector  is  computed  and  this  is  used  to 
determine  h0.  In  most  cases  g  will  be  the  smaller  of  the  norm  of  the  function  and 
the  first  derivative,  but  in  case  one  is  smaller  than  1000*eps  (machine  epsilon)  the 
larger  is  taken,  and  in  no  case  is  g  taken  as  less  than  10'6.  A  6th  order  Runge-Kutta 
method  is  used  by  calling  RK6  to  compute  the  solution  Y4  at  X0+h0/4,  Y2  at 
X0+h0/2  (from  Y4),  and  Y1  at  XO+hO  (from  Y2).  These  are  then  used  to  calculate 
a  Nordsieck  vector  for  computation  of  the  initial  external  stage  vector,  plus  an  extra 
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derivative  for  computation  of  the  appropriate  size  for  the  first  DIMSIM  step.  Both 
the  Nordsieck  vector  and  the  intial  external  stage  vector  are  returned  in  case 
tolerance  is  not  met. 

YIterP:  (real*8  array)  Initial  external  stage  vector,  output  dimensioned  (NEqn,5) 

F,X0,Y0,NEqn:  See  Section  a  above. 

H:  (real  *8)  Output  value  for  calculated  initial  DIMSIM  step  size 

YP:  (real*8  array)  Derivative  of  the  solution  at  X0,  calculated  using  y’=F(X0,Y0),  output 
dimensioned  NEqn 

Tol:  (real*8)  Initial  tolerance  computed  as  described  in  section  a  from  RTol,  Y0  and  ATol, 
input 

Y1,Y1P,Y2,Y2P,Y4:  (real*8  arrays)  Values  computed  using  6th  order  Runge-Kutta 
method  for  solution  at  h0,  h/2,  and  h</4,  and  associated  derivatives  computed  using 
F,  workspace  each  dimensioned  NEqn. 

Y6Den  (real*8  array)  6th  derivative  of  solution  vector  at  X0,  multiplied  times  /$, 
workspace  dimensioned  NEqn. 

YNord:  (real*  8  array)  Nordsieck  vector  at  X0  for  step  size  H,  output  dimensioned 
(NEqn,6) 

RKWork:  See  driver  description  above. 

subroutine  rescale(YIter,H,H01d,FStage,YIterP,NEqn) 

Rescale  takes  a  previous  external  stage  vector  computed  using  a  step  size  HOld  and 
uses  the  DIMSIM  rescaling  algorithm  to  rescale  it  to  become  the  appropriate 
previous  external  stage  vector  for  continuing  the  computation  with  step  size  H. 

YIter:  (real*8  array)  Output  rescaled  external  stage  vector,  dimensioned  (NEqn,5) 

H:  (real*8)  Input  new  step  length 

HOld:  (real*8)  Input  step  length  for  previous  successful  step. 

FStage:  (real*8  array)  Input  vector  of  derivatives  at  stage  points  for  last  step  completed 

YIterP:  (real*  8  array)  Input  previous  external  stage  vector,  calculated  using  step  length 
HOld,  dimensioned  (NEqn, 5) 

NEqn:  See  Section  a  above 

real*8  function  NewStep(H,Tol,ErrEst) 

The  error  estimate  is  compared  with  the  tolerance  to  produce  a  step  change  factor 
modifying  the  step  size  in  preparation  for  the  next  step.  The  local  truncation  error  is 
assumed  to  be  proportional  to  the  6th  power  of  the  step  size.  The  largest  possible 
step  change  is  a  factor  of  2.  A  step  change  safety  factor  equal  to  .75  was  used  for 
DIMEX5  so  that  the  step  size  used  was  .75  of  optimal.  A  larger  factor  of  .9  was 
used  successfully  for  DIMEX2.  Step  size  is  not  changed  if  change  would  be  less 
than  10%  for  DIMEX5  and  5%  for  DIMEX2. 

H:  (real*8)  At  input  contains  old  step  size.  At  output  contains  new  step  size 

Tol:  (real*8)  Input  tolerance  derived  from  RTol,  Y  and  ATol  as  described  in  section  a 
above 

ErrEst:  (real*  8)  Input  current  local  error  estimate 

subroutine  NordInterp(Y,H,Theta,FStage,YIter,NordVec,NEqn) 

The  Nordsieck  vector  NordVec  is  calculated  at  the  current  value  of  XOut,  using  the 
previously  calculated  external  stage  vector  YIter,  the  current  step  size  H,  and  the 
current  vector  of  derivatives  at  stage  points,  FStage.  Theta  then  is  used  to  rescale 
the  Nordsieck  vector  to  correspond  to  the  desired  output  point  and  a  fifthth  order 
Taylor  series  approximation  Y  is  formed. 
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Y:  (real*8  array)  Output  approximation  to  the  solution  at  the  desired  output  point,  length 
NEqn. 

H:  (real*8)  Input  step  size  for  last  successful  integration 

Theta:  (real*  8)  Input  between  0  and  1,  ratio  (X-X01d)/H,  where  X  is  the  desired  output 
point 

FStage:  (real*8  array)  Input  vector  of  derivatives  at  stage  points  for  last  successful 
integration,  dimensioned  (NEqn,5) 

YIter:  (real*8  array)  Input  previous  external  stage  vector,  dimensioned  (NEqn, 5) 

NordVec:  (real*8  array)  Workspace  (output  for  waveform  relaxation)  for  Nordsieck  vector 
at  current  XOut  and  last  step  size  used. 

NEqn;  (integer)  See  above.  Section  a. 

subroutine  RK6(Y,X0,Y0,H,F,NEqn,YP,YP2,YP3,YP4,YP5,YP6,YP7,YTemp) 

A  simple  6th  order  Runge-Kutta  solver,  used  to  provide  solution  values  needed  for 
starting  method.  Solution  Y  for  initial  values  (X0,Y0)  is  found  at  XO+H  for 
problem  Y’=F(X,Y).  The  method  is  described  in  Butcher  (Reference  1 1),  pp.  203- 
205. 

Y:  (real*8  array)  Solution  at  XO+H,  dimensioned  NEqn. 

X0:  (real  *8)  Initial  value  of  independent  variable. 

Y0:  (real*8  array)  Initial  value  of  dependent  variable,  dimensioned  NEqn 

H:  (real*8)  step  size  used 

F:  See  above,  Section  a 

NEqn:  See  above,  Section  a 

YP,YP2,YP3,YP4,YP5,YP6,YP7:  (real*8  arrays)  Workspace  used  for  storing  derivatives 
of  stages,  each  dimensioned  NEqn. 

YTemp:  (real*8  array)  Workspace  used  for  storing  stage  vectors  successfully, 

dimensioned  NEqn 

block  data  matrices 

The  matrices  and  vectors  defining  the  paricular  method  are  set  to  appropriate  values 
using  data  statements  and  stored  in  the  common  region  /Method/. 


TESTING  DIMEX 

The  ultimate  test  of  an  ODE  solver  is  the  accurate  and  efficient  solution  of  a  wide 
variety  of  problems.  The  relatively  easy,  nonstiff  problems  utilized  in  the  well-known  suite 
DETEST  have  provided  a  standard  test  suite  for  ODE  solvers  since  they  were  first 
introduced  in  1972  by  Hull,  et.  al,  (Reference  31).  For  second  order  testing  these  were  run 
using  both  DIMEX2  and  the  well  established  linear  multistep  solver  LSODE  developed  at 
Livermore  by  Hindmarsh  and  his  coworkers  (Reference  32).  LSODE  was  used  in  the 
Adams-Moulton  solver  mode  with  functional  iteration.  This  mode  provides  predictor- 
corrector  solution  of  the  form  P(EC)M  where  M  is  no  greater  than  3  and  is  usually  1 .  This 
provided  the  most  efficient  solution  for  competitive  purposes  on  this  test  suite.  For  fifth 
order  testing  these  were  run  using  both  DIMEX5  and  the  highly  regarded  Runge-Kutta 
code  DOPRI5  (Reference  18).  For  the  easy  problems  of  DETEST,  LSODE  is  very 
efficient  as  measured  by  steps  taken  and  function  evaluations,  but  the  sphere  of 
applicability  for  DIMEX5  would  be  similar  to  that  of  DOPRI5  and  it  is  hoped  that  a 
comparison  on  simpler  problems  would  provide  some  indication  of  relative  performance  on 
the  harder  problems  that  do  not  find  their  way  into  standard  benchmark  suites.  The 
problems  are  divided  into  classes  A-E,  with  Class  A  consisting  of  single  equations,  Class 
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B  of  small  systems  of  2  or  3  equations.  Class  C  of  moderate  systems  of  10  to  51 
equations,  Class  D  of  orbit  equations,  and  Class  E  of  higher  order  (second  order) 
equations.  Application  areas  represented  include  the  simple  negative  exponential,  a  simple 
Riccati  equation,  an  oscillatory  problem,  logistics  and  spiral  curves,  conflicting  population 
growths,  chemical  reactions,  rigid  body  motion,  radioactive  decay  chains,  heat  equation,  5 
body  problem,  orbital  motion,  Bessel’s  equation,  Van  der  Pol’s  equation,  Duffing’s 
equation,  falling  body  and  linear  pursuit.  All  are  integrated  over  an  interval  from  0  to  20. 
Only  a  few  are  provided  with  exact  solutions,  and  accuracy  testing  was  conducted  primarily 
using  significantly  tighter  tolerances  with  the  same  solver.  Rather  than  adapting  DIMEX, 
DOPRI5  and  LSODE  to  fit  with  the  DETEST  software,  the  functions  were  simply  coded 
and  used  with  a  specially  designed  calling  routine  that  provided  the  desired  output 
information.  Although  number  of  steps  is  important,  even  more  important  is  the  number  of 
function  evaluations,  which  is  a  good  measure  of  the  amount  of  work  required.  A  timing 
of  the  solution  process  would  also  give  a  measure  of  the  overhead,  which  has  been  noted  to 
be  significant  for  LSODE  (Reference  18).  But  the  DIMEX  codes  have  not  been  carefully 
optimized  yet  and  so  this  comparison  was  not  carried  out. 

Some  aspects  of  the  comparison  process  utilized  should  be  noted.  First,  LSODE  uses 
an  Adams-Moulton  solver  and  can  only  be  set  to  restrict  to  a  specified  order  up  to  12  that 
will  not  be  exceeded.  Thus  for  second  order  testing  it  is  possible  that  some  steps  were  first 
order.  The  order  is  chosen  to  minimize  error  and  hence  maximizes  step  length. 
Furthermore,  in  conducting  these  tests  it  was  observed  that  LSODE  applies  relative 
tolerance  to  each  separate  solution  component,  while  DIMEX  applies  relative  tolerance  to 
the  norm  of  the  solution.  Thus  the  only  relative  tolerance  comparisons  are  for  tests  A1-A5, 
since  only  a  single  equation  is  used.  Otherwise,  absolute  tolerance  was  used  for  both. 
Finally,  exact  solutions  were  used  for  comparison  purposes  only  for  A 1-4,  while  high 
tolerances  of  10'13  to  10'15  were  used  to  produce  a  “true”  value  of  the  local  solution  and  the 
end  point  solution. 

For  fifth  order  testing  was  conducted  only  for  the  more  appropriate  higher  tolerances 
10'6,  10'9,  and  10'12.  On  the  other  hand,  second  order  methods  often  take  a  very  long  time 
to  integrate  at  tight  tolerances  and  so  10'12  was  not  used  for  second  order,  while  the  loose 
tolerance  of  10"3  was  an  appropriate  condition  for  testing  a  second  order  solver. 

A  summary  of  results  is  shown  in  the  Tables  22  through  71.  Note  that  an  entry  of  “**” 
indicates  that  a  measurement  was  not  made.  It  was  not  convenient  to  obtain  data  from 
LSODE  on  tolerance  misses  and  times  deceived,  and  it  did  not  seem  to  be  important  to 
obtain  that  information.  DOPRI5  provided  number  of  tolerance  misses  but  it  was  not 
convenient  to  obtain  the  number  of  times  deceived.  And  DIMEX5  did  not  conveniently 
provide  a  sufficiently  accurate  solution  at  tolerance  10'15  to  determine  the  accuracy  of 
individual  steps  as  required  in  computing  number  of  steps  deceived  at  tolerance  10"!2. 
Also,  “AM”  designates  the  Adams-Moulton  integrator  of  LSODE. 
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DIMEX2  Test  Results 


TABLE  22.  Order  2  Problem  A1  Test. 


Tolerance 

10' 

9 - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

121 

107 

1222 

1049 

12231 

#  func  calls 

126 

116 

1227 

1055 

12235 

#  tol  misses 

1 

** 

1 

** 

0 

** 

#  deceived 

0 

0 

0 

** 

0 

** 

end  relerr 

7.9e-2 

1.3e-l 

8.8e-4 

6.8e-4 

8.9e-6 

6.1e-6 

TABLE  23.  Order  2  Problem  A2  Test. 


Tolerance 

10 

3 - 

5 - 

io-> 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

27 

29 

239 

217 

2356 

#  func  calls 

31 

32 

243 

220 

2360 

#  tol  misses 

0 

** 

0 

** 

0 

** 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

3.4e-3 

4.3e-3 

4.8e-5 

3.5e-5 

5.0e-7 

3.3e-7 

TABLE  24.  Order  2  Problem  A3  Test. 


Tolerance 

■H 

5 - 

np 

method 

DM2 

AME2 

DM2 

AM2 

DM2 

AM2 

#  steps 

131 

108 

1160 

950 

11485 

9430 

#  func  calls 

171 

132 

1199 

1002 

11521 

9523 

#  tol  misses 

33 

** 

31 

** 

32 

** 

#  deceived 

2 

** 

1 

** 

0 

** 

end  relerr 

1.3e-3 

4.6e-3 

4.2e-5 

4.6e-5 

5.8e-7 

4.2e-7 

TABLE  25.  Order  2  Problem  A4  Test. 


Tolerance 

"3 - 1 

ms 

5 - 

lO" 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

19 

18 

169 

144 

1650 

1377 

#  func  calls 

25 

26 

178 

156 

1659 

1386 

#  tol  misses 

2 

** 

5 

** 

5 

** 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

2.5e-4 

4.1e-4 

1.6e-6 

1.9e-7 

7.0e-9 

5.8e-9 
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TABLE  26.  Order  2  Problem  A5  Test. 


Tolerance 

10' 

3  " 

5 

10 

-y 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

24 

21 

185 

162 

1801 

1539 

#  func  calls 

37 

33 

197 

179 

1815 

1566 

#  tol  misses 

9 

** 

8 

** 

10 

** 

#  deceived 

0 

** 

11 

** 

9 

** 

end  relerr 

1.5e-l 

2.7e-2 

1.6e-3 

8.4e-4 

1.6e-5 

1.0e-5 

TABLE  27.  Order  2  Problem  B 1  Test. 


Tolerance 

■J  "" 

HHHBI 

b 

10 

=9 - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

244 

191 

2194 

1724 

21944 

17060 

#  func  calls 

290 

224 

2199 

1766 

21949 

17102 

#  tol  misses 

42 

** 

1 

** 

1 

** 

#  deceived 

32 

** 

1 

** 

0 

** 

end  relerr 

4.9e-3 

5.4e-l 

8.6e-5 

5.4e-4 

l.le-6 

3.5e-6 

TABLE  28.  Order  2  Problem  B2  Test. 


Tolerance 

b 

10 

-9 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

92 

90 

275 

239 

2423 

#  func  calls 

101 

151 

284 

276 

2427 

1936 

#  tol  misses 

5 

** 

5 

** 

0 

** 

#  deceived 

0 

** 

3 

** 

0 

** 

end  relerr 

1.2e-4 

4.2e-5 

2.4e-7 

2.6e-8 

2.5e-10 

l.le-9 

TABLE  29.  Order  2  Problem  B  3  Test. 


Tolerance 

S - 

10 

-y 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

|  AM2 

#  steps 

37 

45 

259 

209 

2548 

#  func  calls 

45 

72 

263 

212 

2552 

#  tol  misses 

4 

** 

0 

** 

0 

** 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

4.5e-4 

3.7e-4 

1.7e-5 

1.7e-5 

1.9e-7 

1.4e-7 
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TABLE  30.  Order  2  Problem  B4  Test. 


Tolerance 

b 

10 

r9 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

191 

143 

1882 

1356 

18833 

#  func  calls 

197 

150 

1886 

1363 

18837 

13713 

#  tol  misses 

2 

** 

0 

** 

0 

** 

#  deceived 

8 

** 

0 

** 

0 

** 

end  relerr 

1.9e-l 

1.0e-l 

1.8e-3 

1.7e-3 

1.8e-5 

1.6e-5 

TABLE  31.  Order  2  Problem  B5  Test. 


Tolerance 

•3  " . 

io- 

b 

10 

^ - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

114 

87 

1048 

788 

10378 

7686 

#  func  calls 

133 

150 

1066 

821 

10384 

7724 

#  tol  misses 

15 

** 

14 

** 

3 

** 

#  deceived 

10 

** 

0 

** 

0 

** 

end  relerr 

1.6e-2 

l.le-1 

3.7e-4 

6.1e-4 

3.1e-6 

5.0e-6 

TABLE  32.  Order  2  Problem  Cl  Test. 


Tolerance 

10‘ 

■5 

5 - 

10 

^ - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

54 

50 

356 

247 

3523 

2383 

#  fimc  calls 

60 

80 

360 

250 

3527 

2386 

#  tol  misses 

2 

** 

0 

** 

0 

** 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

2.1e-4 

2.1e-4 

1.3e-5 

2.6e-5 

1.4e-7 

1.8e-7 

TABLE  33.  Order  2  Problem  C2  Test. 


Tolerance 

•5 

10" 

5 - 

10 

-y 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

268 

225 

394 

325 

2961 

1949 

#  func  calls 

276 

406 

403 

439 

2970 

1995 

#  tol  misses 

4 

** 

5 

** 

5 

** 

#  deceived 

0 

** 

2 

** 

0 

** 

end  relerr 

1.5e-4 

2.7e-5 

3.9e-7 

1.3e-7 

3.5e-10 

6.9e-10 
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TABLE  34.  Order  2  Problem  C3  Test. 


Tolerance 

IHHHQ9 

■5 

5 - 

10 

-y 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

116 

112 

259 

202 

2352 

1567 

#  func  calls 

123 

205 

266 

245 

2356 

1570 

#  tol  misses 

3 

** 

3 

** 

0 

** 

#  deceived 

1 

** 

0 

** 

0 

** 

end  relerr 

5.5e-2 

9.5e-3 

1.9e-4 

1.2e-4 

4.5e-6 

5.0e-6 

TABLE  35.  Order  2  Problem  C4  Test. 


Tolerance 

MHHI9 

1 - 

b 

10 

-y 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

117 

93 

259 

165 

2355 

1199 

#  func  calls 

124 

167 

266 

199 

2359 

#  tol  misses 

3 

** 

3 

** 

0 

** 

#  deceived 

1 

** 

3 

** 

1 

** 

end  relerr 

1.4e-2 

2.5e-2 

1.8e-4 

2.0e-4 

4.1e-6 

7.0e-6 

TABLE  36.  Order  2  Problem  C5  Test. 


Tolerance 

■3 - 

3 

10 

-y - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

30 

16 

287 

157 

2865 

1536 

#  func  calls 

34 

21 

291 

159 

2869 

1538 

#  tol  misses 

0 

** 

0 

** 

0 

** 

#  deceived 

0 

** 

0 

** 

1 

** 

end  relerr 

1.4e-3 

4.9e-3 

1.2e-5 

1.2e-5 

l.le-7 

1.9e-7 

TABLE  37.  Order  2  Problem  D1  Test. 


Tolerance 

HHflDSI 

5 

Bun 

S' 

10 

-9 - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

130 

115 

1255 

937 

12541 

9418 

#  func  calls 

135 

168 

1259 

944 

12545 

9426 

#  tol  misses 

1 

** 

0 

** 

1 

** 

#  deceived 

6 

** 

0 

** 

0 

** 

end  relerr 

7.6e-l 

2.1e+0 

4.0e-3 

4.6e-3 

3.4e-5 

9.4e-6 
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TABLE  38.  Order  2  Problem  D2  Test. 


Tolerance 

10 

3 - 

5 - 

10 

r5 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

155 

129 

1483 

1071 

14816 

10659 

#  func  calls 

165 

199 

1487 

1091 

14820 

10680 

#  tol  misses 

6 

** 

0 

** 

0 

** 

#  deceived 

11 

** 

0 

** 

0 

** 

end  relerr 

6.3e-l 

2.0e+0 

5.7e-3 

2.5e-3 

5.8e-5 

2.1e-5 

TABLE  39.  Order  2  Problem  D3  Test. 


Tolerance 

10‘ 

3 - 

Mil 

5 - 

10 

^ - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

202 

152 

1799 

1289 

17937 

12746 

#  func  calls 

236 

230 

1803 

1322 

17941 

12779 

#  tol  misses 

30 

** 

0 

** 

0 

** 

#  deceived 

6 

** 

0 

** 

0 

** 

end  relerr 

7.1e-l 

2.2e+0 

7.6e-3 

1.7e-3 

8.2e-5 

4.1e-5 

TABLE  40.  Order  2  Problem  D4  Test. 


Tolerance 

10‘ 

•4 

5 - 

10 

- 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

269 

188 

2294 

1634 

22896 

16255 

#  func  calls 

330 

256 

2298 

1683 

22900 

16306 

#  tol  misses 

57 

** 

0 

** 

0 

** 

#  deceived 

8 

** 

0 

** 

0 

m* 

end  relerr 

7.0e-l 

1.8e+0 

8.9e-3 

3.0e-3 

9.3e-5 

5.6e-5 

TABLE  41.  Order  2  Problem  D5  Test. 


Tolerance 

10 

•4 . . 

5 

10 

- 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

414 

280 

3457 

2464 

34502 

24212 

#  func  calls 

521 

363 

3461 

2545 

34506 

24296 

#  tol  misses 

103 

** 

0 

** 

0 

** 

#  deceived 

11 

** 

0 

** 

1 

** 

end  relerr 

5.4e-l 

1.2e+0 

9.3e-3 

5.2e-3 

9.3e-5 

7.4e-5 
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TABLE  42.  Order  2  Problem  El  Test. 


Tolerance 

•i 

5 - 

10 

^ - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

74 

62 

738 

611 

7372 

6021 

#  func  calls 

78 

117 

742 

614 

7376 

6024 

#  tol  misses 

0 

** 

0 

** 

0 

** 

#  deceived 

8 

** 

0 

** 

0 

** 

end  relerr 

2.7e-l 

2.6e-l 

2.0e-3 

1.8e-3 

2.0e-5 

1.5e-5 

TABLE  43.  Order  2  Problem  E2  Test. 


Tolerance 

•5 

5 

10 

^ - 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

282 

219 

2649 

20  io 

26490 

#  func  calls 

318 

264 

2657 

2061 

26495 

#  tol  misses 

32 

** 

4 

** 

1 

** 

#  deceived 

28 

** 

0 

** 

0 

** 

end  relerr 

2.2e-2 

1.3e-2 

1.9e-4 

1.5e-4 

1.9e-6 

1.6e-6 

TABLE  44.  Order  2  Problem  E3  Test. 


Tolerance 

* - 

10" 

5 

10 

-y - 

method 

DM2 

AM2 

DIM2 

AM2 

DM2 

AlvG 

#  steps 

281 

213 

2714 

2064 

27128 

#  func  calls 

304 

230 

2718 

2086 

27133 

#  tol  misses 

19 

** 

0 

** 

I 

** 

#  deceived 

3 

** 

0 

** 

0 

** 

end  relerr 

9.7e-2 

1 ,4e- 1 

1.2e-3 

1.2e-3 

l.le-5 

1.3e-5 

TABLE  45.  Order  2  Problem  E4  Test. 


5olerance 

•3 

IHKI 

10 

-y . 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

i  AM2 

#  steps 

16 

14 

138 

104 

1365 

#  func  calls 

22 

23 

142 

109 

1369 

#  tol  misses 

2 

** 

0 

** 

0 

** 

#  deceived 

1 

** 

5 

** 

0 

** 

end  relerr 

2.8e-4 

3.7e-4 

4.3e-6 

3.7e-6 

4.3e-8 

3.8e-8 
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TABLE  46.  Order  2  Problem  E5  Test. 


Tolerance 

3 - 

5 - 

- 

method 

DM2 

AM2 

DM2 

AM2 

DM2 

AM2 

#  steps 

31 

28 

280 

209 

2804 

wmssnm 

#  func  calls 

32 

52 

285 

220 

2809 

#  tol  misses 

1 

** 

1 

** 

1 

** 

#  deceived 

23 

** 

0 

** 

0 

** 

end  relerr 

4.9e-3 

3.5e-3 

4.7e-5 

4.2e-5 

4.7e-7 

4.4e-7 

DIMEX5  TEST  RESULTS 


TABLE  47.  Order  5  Problem  A1  Test. 


Tolerance 

•0 

9 - 

10 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

97 

94 

288 

361 

882 

1426 

#  func  calls 

531 

566 

1496 

2168 

4436 

8558 

#  tol  misses 

5 

0 

7 

0 

1 

0 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

l.le-5  1 

3.9e-6 

5.0e-8 

3.3e-9 

1.7e-10 

3.1e-12 

TABLE  48.  Order  5  Problem  A2  Test. 


Tolerance 

■S - 

■am 

5 - 

Hr7 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

38 

22 

104 

67 

320 

246 

#  func  calls 

226 

134 

546 

404 

1621 

1478 

#  tol  misses 

3 

0 

1 

b 

0 

0 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

8.7e-7 

4.4e-7 

5.1e-9 

2.8e-10 

1.7e-ll 

1.9e-13 

TABLE  49.  Order  5  Problem  A3  Test. 


Tolerance 

9 - ! 

JqTTT 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

■MM 

136 

88 

426 

286 

1348 

1070 

#  func  calls 

721 

530 

2166 

1718 

6776 

6422 

#  tol  misses 

4 

15 

3 

19 

3 

19 

#  deceived 

o 

** 

0 

** 

b 

** 

end  relerr 

4.3e-6 

2.1e-6 

8.8e-9 

1.7e-9 

2.6e-l  1 

2.1e-12 

102 


NAWCWPNS  TP  8340 


TABLE  50.  Order  5  Problem  A4  Test. 


Tolerance 

■HHB9 

0 

10': 

? - 

lO'12 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

23 

18 

64 

56 

197 

205 

#  func  calls 

146 

110 

346 

338 

1016 

1232 

#  tol  misses 

2 

1 

1 

2 

2 

1 

#  deceived 

0 

** 

0 

** 

0 

** 

end  relerr 

1.9e-7 

7.2e-8 

3.2e-10 

5.1e-ll 

9.8e-14 

6.9e-13 

TABLE  5 1 .  Order  5  Problem  A5  Test. 


Tolerance 

O 

- W - 

- W2 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

28 

15 

78 

51 

241 

190 

#  func  calls 

171 

92 

411 

308 

1226 

1142 

#  tol  misses 

2 

6 

0 

1 

0 

1 

#  deceived 

2 

** 

1 

** 

** 

** 

end  relerr 

5.4e-5 

6.5e-6 

2.3e-7 

7.6e-9 

3.5e-10 

8.9e-12 

TABLE  52.  Order  5  Problem  B 1  Test. 


Tolerance 

S - 

10' 

9 

10-ii 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

233 

165 

683 

583 

2147 

2311 

#  func  calls 

1201 

992 

3436 

3500 

10756 

13868 

#  tol  misses 

3 

i7 

0 

3 

0 

0 

#  deceived 

2 

** 

1 

** 

** 

** 

end  relerr 

2.1e-4 

1.9e-5 

9.9e-7 

2.5e-9 

3.3e-9 

2.3e-12 

TABLE  53.  Order  5  Problem  B2  Test. 


Tol  (abs) 

HHUS 

0 

10' 

? - 

lO1 

2 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

87 

43 

134 

123 

354 

450 

#  func  calls 

486 

260 

711 

740 

1806 

2702 

#  tol  misses 

6 

2 

4 

0 

3 

0 

#  deceived 

6 

** 

0 

** 

** 

** 

end  relerr 

3.6e-7 

4.7e-7 

2.9e-10 

6.3e-ll 

2.1e-13 

9.7e-14 
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TABLE  54.  Order  5  Problem  B3  Test. 


Tol  (abs) 

10 

■b 

hh 

9 - 

10"' 

method 

DLM5 

DOPRI5 

DIM5 

DOPRI5 

D1M5 

DOPRI5 

#  steps 

49 

33 

112 

108 

344 

400 

#  func  calls 

266 

200 

596 

650 

1741 

2402 

#  tol  misses 

4 

0 

3 

0 

0 

** 

#  deceived 

i 

** 

0 

** 

** 

** 

end  relerr 

2.7e-7 

9.7e-8 

2.1e-9 

1.9e-10 

7.4e-12 

1.4e-13 

TABLE  55.  Order  5  Problem  B4  Test. 


Tol  (abs) 

■t> 

9 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

159 

108 

494 

424 

1561 

1682 

#  func  calls 

816 

650 

2491 

2546 

7826 

10094 

#  tol  misses 

0 

0 

0 

0 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

3.1e-6 

3.4e-6 

5.8e-9 

8.1e-10 

2.8e-ll 

8.9e-13 

TABLE  56.  Order  5  Problem  B5  Test. 


Tol  (abs) 

10 

-5 - 

9 - 

10"' 

method 

DM5 

DOPRI5 

DM5 

DOPR15 

DM5 

DOPRI5 

#  steps 

137 

82 

422 

284 

1327 

1123 

#  func  calls 

706 

494 

2136 

1706 

6661 

6740 

#  tol  misses 

0 

7 

1 

1 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

6.8e-6 

3.0e-6 

3.2e-8 

7.7e-9 

8.9e-l  1 

9.0e-12 

TABLE  57. 

Order  5  Problem  Cl  Test. 

Tol  (abs) 

b 

io-; 

9 

10"' 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

51 

39 

132 

134 

411 

508 

#  func  calls 

281 

236 

681 

806 

2076 

3050 

#  tol  misses 

1 

0 

0 

0 

0 

0 

#  deceived 

9 

** 

9 

** 

** 

** 

end  relerr 

l.le-7 

1.7e-7 

1.5e-9 

4.8e-10 

5.7e-12 

5.5e-13 
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TABLE  58.  Order  5  Problem  C2  Test. 


Tol  (abs) 

5 - 

10- 

1 

lO^ 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

198 

68 

214 

142 

466 

496 

#  func  calls 

1086 

410 

1126 

854 

2381 

2978 

#  tol  misses 

15 

2 

7 

2 

6 

0 

#  deceived 

17 

** 

43 

** 

** 

** 

end  relerr 

2.6e-8 

7.2e-7 

7.8e~ll 

2.2e-10 

2.2e-13 

2.7e-13 

TABLE  59.  Order  5  Problem  C3  Test. 


Tol  (abs) 

Z - 

- To™ 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

98 

41 

122 

108 

334 

398 

#  func  calls 

536 

248 

646 

650 

1691 

2390 

#  tol  misses 

5 

i 

3 

0 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

1.4e-5 

3.2e-5 

2.2e-8 

1.4e-8 

1.9e-10 

1.8e-ll 

TABLE  60.  Order  5  Problem  C4  Test. 


Tol  (abs) 

z - 

3 - 

To™ 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

104 

39 

132 

94 

334 

340 

#  func  calls 

561 

236 

696 

566 

1691 

2042 

#  tol  misses 

4 

1 

3 

0 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

1.0t-l 

3.5e-5 

7.6e-9 

3.0e-8 

1.5e-10 

3.4e-ll 

TABLE  61.  Order  5  Problem  C5  Test. 


Tolerance 

10' 

■ 

r™“ 

10-  iz 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

14 

15 

32 

55 

94 

211 

#  func  calls 

91 

92 

181 

332 

491 

1268 

#  tol  misses 

0 

0 

0 

2 

0 

4 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

1.8e-6 

1.8e-6 

6.7e-9 

8.3e-10 

1.8e-ll 

5.7e-13 
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TABLE  62.  Order  5  Problem  D1  Test. 


Tol  (abs) 

10' 

3 - 

9 - ; 

1012 

method 

DIM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

119 

83 

373 

323 

1177 

1276 

#  func  calls 

616 

500 

1891 

1940 

5911 

7658 

#  tol  misses 

0 

0 

1 

6 

1 

0 

#  deceived 

0 

** 

1 

** 

** 

** 

end  releir 

1.3e-4 

3.1e-4 

8.4e-7 

1.3e-7 

2.5e-9 

1.7e-10 

TABLE  63.  Order  5  Problem  D2  Test. 


Tol  (abs) 

■O 

9 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

163 

98 

512 

350 

1617 

1382 

#  func  calls 

836 

590 

2581 

2102 

8106 

8294 

#  tol  misses 

6 

6 

0 

0 

6 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

1.5e-4 

4.3e-4 

4.5e-7 

1.2e-7 

1.4e-9 

1.5e-9 

TABLE  64.  Order  5  Problem  D3  Test. 


Tol  (abs) 

■IHH9 

^ - : 

BflHSl 

9 

10[ 

7 - 

method 

DM5 

DOPRI5 

|  DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

212 

131 

671 

414 

2121 

1638 

#  func  calls 

1081 

788 

3376 

2486 

10631 

9830 

#  tol  misses 

0 

20 

0 

0 

1 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

4.0e-4 

2.6e-4 

1.0e-6 

9.1e-8 

2.8e-9 

l.le-10 

TABLE  65.  Order  5  Problem  D4  Test. 


Tol  (abs) 

■HHB3 

■6  " 1 

10T 

7 - 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

281 

174 

889 

514 

2807 

2035 

#  func  calls 

1426 

1046 

4471 

3086 

14061 

12212 

#  tol  misses 

0 

34 

1 

0 

1 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

6.8e-4 

2.0e-4 

1.5e-6 

8.0e-8 

4.3e-9 

9.7e-ll 
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TABLE  66.  Order  5  Problem  D5  Test. 


Tol  (abs) 

•6 

9 - 

W2 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

422 

256 

1331 

728 

4201 

2884 

#  func  calls 

2136 

1538 

6686 

4370 

21056 

17306 

#  tol  misses 

1 

54 

0 

0 

0 

0 

#  deceived 

1 

** 

1 

** 

** 

** 

end  relerr 

67.1-5 

1.3e-4 

9.1e-8 

6.0e-8 

2.2e-9 

6.7e-ll 

TABLE  67.  Order  5  Problem  El  Test. 


Tolerance 

S - 

9 

10'1Z 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

94 

112 

286 

426 

900 

1683 

#  func  calls 

491 

674 

1451 

2558 

4521 

10100 

#  tol  misses 

0 

10 

0 

11 

0 

12 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

l.le-5 

2.2e-6 

4.3e-8 

2.7e-9 

1.3e-10 

2.8e-12 

TABLE  68.  Order  5  Problem  E2  Test. 


Tol  (abs) 

9 

10'iz 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

303 

183 

916 

637 

2922 

2503 

#  func  calls 

1551 

1100 

4601 

3824 

14631 

15020 

#  tol  misses 

3 

18 

0 

6 

0 

5 

#  deceived 

2 

** 

1 

** 

** 

** 

end  relerr 

5.0e-6 

2.0e-4 

4.3e-9 

6.1e-8 

7.3e-12 

2.6e-ll 

TABLE  69.  Order  5  Problem  E3  Test. 


Tol  (abs) 

Z - 

101 

2 

method 

DM5 

DOPRI5 

DM5 

DOPRI5 

DM5 

DOPRI5 

#  steps 

220 

137 

695 

489 

2198 

1916 

#  func  calls 

1121 

824 

3501 

2936 

11011 

11498 

#  tol  misses 

0 

10 

1 

4 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

8.8e-6 

1.6e-4 

4.4e-8 

4.9e-8 

4.0e-10 

2.3e-ll 
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TABLE  70.  Order  5  Problem  E4  Test. 


Tol  (abs) 

10 

■0 

io- 

y 

- W1 

method 

DIM5 

DOPRI5 

DIM5 

DOPRI5 

DIM5 

DOPRI5 

#  steps 

19 

17 

47 

46 

145 

167 

#  func  calls 

121 

104 

256 

278 

746 

1004 

#  tol  misses 

1 

1 

0 

0 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

3.9e-8 

1.8e-7 

1.3e-10 

1.5e-10 

1.3e-13 

1.7e-13 

TABLE  71.  Order  5  Problem  E5  Test. 


Tol  (abs) 

■t> 

9 - 

- 

method 

DIM5 

DOPRI5 

DIM5 

DOPRI5 

DIM5 

DOPRI5 

#  steps 

25 

24 

78 

61 

242 

233" 

#  func  calls 

146 

146 

411 

368 

1231 

1400 

#  tol  misses 

0 

5 

0 

0 

0 

0 

#  deceived 

0 

** 

0 

** 

** 

** 

end  relerr 

1.3e-6 

1.6e-8 

4.3e-9 

1.5e-ll 

1.2e-ll 

1.2e-14 

DISCUSSION  OF  TEST  RESULTS 

In  scanning  the  results  for  second  order  testing,  it  is  evident  that  the  number  of  steps  is 
usually  10  to  20%  higher  for  the  DIMEX2  integrator  in  comparison  with  LSODE.  But 
since  LSODE  has  been  extensively  refined  over  a  period  of  nearly  20  years,  this  is  not  bad 
for  a  newly  developed  research  code.  A  different  DIMSIM  with  a  smaller  error  constant 
would  produce  fewer  steps  with  otherwise  the  same  implementation.  On  the  other  hand, 
for  these  relatively  easy  problems,  LSODE  minimizes  the  number  of  function  iterations  to 
the  extent  that  often  the  average  number  of  function  evaluations  required  is  just  over  1  per 
step.  DIMEX2  significantly  outperforms  LSODE  at  loose  tolerances  on  problems  A2,  A4, 
B2,  B3,  B5,  Cl,  C2,  C3,C4,  El,  E4,  and  E5.  The  orbital  problems  D1  and  D2  might  also 
be  included,  but  the  accuracy  of  both  solvers  was  inadequate  to  make  these  true  tests.  The 
operation  of  DIMEX2  shows  that  although  it  is  not  really  a  competitive  method  for  simple 
problems  with  tighter  tolerances,  the  implementation  is  essentially  correct  and  appropriate. 
And  it  points  to  the  possibility  of  developing  a  truly  competitive  DIMSIM  second  order 
explicit  solver,  since  the  number  of  function  evaluations  per  step  is  consistently  1  even  at 
looser  tolerances. 

A  calculation  of  the  stability  region  for  the  FASAL  implementation  was  performed. 
The  following  graph  shows  that  the  portion  of  the  negative  real  axis  included  was  reduced 
by  1/3.  However,  it  was  found  that  method  parameter  could  be  chosen  to  find  another 
Type  1  second  order  DIMSIM  with  a  more  favorable  stability  region.  The  longest  portion 
of  the  negative  real  axis  that  could  be  obtained  is  also  shown  and  is  unreduced  from  the 
original  DIMSIM,  although  the  portion  off  the  real  axis  is  seriously  reduced  for  this  case  in 
Figure  14.  It  is  felt  that  more  can  be  done  to  produce  desirable  stability  regions  for  higher 
orders  with  more  parameters. 
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FIGURE  14.  Comparison  of  Stability  Regions  for  Second  Order. 


The  fifth  order  results  show  very  competitive  results.  DIMEX5  often  beats  DOPRI5, 
and  otherwise  provides  roughly  comparable  results,  and  this  despite  the  larger  error 
constant  (1/720  versus  1/3600,  a  factor  of  5)  and  a  slightly  smaller  stability  region.  This  is 
directly  attributable  to  the  savings  of  one  function  evaluation  per  step  with  the  DIMSIM. 
But  this  does  not  show  the  complete  picture,  because  dense  output  is  typically  required  in 
real  applications,  not  just  a  race  to  the  end  of  the  interval.  And  the  DIMSIM  requires  no 
additional  function  evaluations  or  system  solves  to  provide  interpolated  output  at  5th  order, 
while  DOPRI5  requires  two  additional  stages  to  provide  5th  order  dense  output  (Reference 
18).  Thus  when  dense  output  requirements  are  considered,  DIMEX5  in  reality 
significantly  outperforms  DOPRI5  for  much  of  this  test  set  and  is  rarely  beaten.  And 
improved  DIMSIMs  may  be  found  with  smaller  error  constants  and  larger  stability  regions. 
Finally,  the  difference  at  higher  orders  is  expected  to  be  much  greater.  For  example,  at  8th 
order,  the  Runge-Kutta-based  DOP853  code  provides  an  8th  order  method  with  an 
interpolant  of  7th  order  with  an  effective  number  of  15  function  evaluations  per  step 
(Reference  18),  while  an  8th  order  DIMSIM  requires  only  8  function  evaluations  per  step 
and  provides  an  interpolant  of  8th  order. 

A  note  should  be  provided  concerning  FASAL  operation.  The  reported  tests  were  first 
conducted  using  the  FASAL  property  only  when  step  size  had  not  changed  since  the  last 
step,  and  for  most  problems  the  number  of  function  evaluations  was  reduced  by  the 
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expected  20%.  An  attempt  was  made  to  allow  this  for  all  steps  after  the  first,  and  in  most 
problems  the  results  were  satisfactory  and  the  number  of  function  evaluations  significantly 
reduced.  However,  for  perhaps  5  or  6  of  the  problems  the  quality  of  the  error  estimate 
deteriorated  significantly,  resulting  in  unacceptably  large  numbers  of  times  the  error 
estimator  was  deceived,  reducing  accuracy,  and  causing  up  to  one  third  more  function 
evaluations.  Further  examination  showed  that  for  these  problems,  not  using  FASAL 
produced  the  best  results.  Table  72  shows  what  happened  for  Problem  C2  with  tolerance 
10'9. 


TABLE  72.  Order  5  With  and  Without  FASAL  for  Problem  C2  Test. 


Use  of  FASAL 

None 

Constant  Step 

All  After  First 

#Steps 

234 

357 

^  590 

#func  calls 

1286 

1633 

2386 

#tol  misses 

19 

7 

1 

#deceived 

81 

209 

430 

end  relerr 

3.0e-10 

2.5e-ll 

2.0e-10 

It  was  originally  thought  that  this  is  due  to  larger  higher  order  terms  that  are  present  in 
these  problems,  and  this  may  be  an  important  factor.  However,  a  calculation  showed  that 
the  stability  region  for  order  5  became  extremely  limited  with  FASAL  implementation  (see 
Figure  15). 
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FIGURE  15.  Order  5  FASAL  Implementation  Stability  Region. 


In  Figure  15  z  values  are  shown  in  a  neighborhood  of  the  origin  where  the  eigenvalue  of 
largest  modulus  w  of  the  stability  matrix  has  Iwl  =  1.  The  dimensions  of  this  region  are 
smaller  by  an  order  of  magnitude  as  compared  with  the  original  DIMSIM.  It  is  evident  that 
stability  for  the  FASAL  implementation  must  be  considered  when  deriving  DIMSIMs  if  this 
implementation  efficiency  is  to  be  utilized!  As  a  result  of  the  small  stability  region,  a 
FASAL  implementation  was  not  used  in  the  testing  reported  above. 

For  the  DIMSIM  itself,  the  error  constant  for  the  method  is  rather  small,  and  at  times 
the  sixth  order  term  might  be  dominated  by  the  seventh  order  terms.  If  this  is  true,  a 
different  method  might  not  have  this  problem.  More  study  of  the  estimation  of  higher  order 
terms  and  of  the  use  of  FASAL  is  needed. 


FUTURE  DEVELOPMENT 

It  might  be  noted  that  DIMEX5  frequently  generates  exactly  1  tolerance  miss.  This  is 
usually  generated  on  the  second  step,  when  the  step  size  calculated  for  the  first  step  is  used 
for  the  second  step,  which  has  a  larger  associated  error.  A  future  version  of  the  code  will 
incorporate  an  appropriate  modification  similar  to  the  approach  used  in  DIMEX2.  Other 
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improvements  and  optimizations  are  anticipated  for  later  development  and  suggestions  and 
bug  reports  are  requested.  It  should  also  be  noted  that  DIMEX5  uses  one  example  of  a 
fifth  order  Type  1  DIMSIM.  It  is  expected  that  in  the  future  other  methods  of  this  class  will 
be  found  with  different  error  constants  and  relatively  smaller  coefficients  for  higher  order 
terms,  and  these  may  enable  significant  performance  improvements.  Stability  for  FASAL 
implementations  will  also  be  utilized  in  deriving  future  methods.  And  although  Runge- 
Kutta  stability  regions  are  a  design  choice  here,  they  are  not  necessarily  essential  and 
stability  regions  even  larger  than  those  of  Runge-Kutta  methods  may  be  found.  And  the 
family,  which  now  includes  only  orders  2  and  5,  may  already  be  extended  to  include 
methods  through  order  8.  It  is  hoped  that  in  time  a  single  variable  order  code  may  be 
developed  incorporating  the  entire  family  in  an  optimal  way. 
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Appendix 

DIMEX5  Code 


subroutine  dim5 1  (X0, Y0, X,Y,F,NEqn,XF,H,ATol,RTol, Starting, XOut, 

1  NCalls,NMisses,Diag,ErrEst,Work) 

£**********:{:******;*::{:**************************************************** 

**** 

c 

c  dim51  calls  driver  which  controls  the  integration  process.  Main 
c  function  of  dim5 1  is  to  divide  up  work  array  so  its  use  is  transparent 

c  to  user.  Work  array  is  dimensioned  at  least  52*NEqn+2. 

c  Parameter  definitions: 

c  X0:  (real*  8)  Initial  value  of  independent  variable  at  start  of  integration, 
c  input 

c  Y0:  (real*8  array)  Initial  value  of  solution  at  start  of  integration, 
c  dimensioned  NEqn,  input 

c  X:  (real*8)  Desired  output  point  if  Diag  is  .false.,  otherwise  set  to  same 
c  as  XF,  input 

c  Y:  (real*8)  Computed  solution  at  X  if  Diag  is  .false.,  otherwise  at  XOut, 
c  output  only 

c  F:  (subroutine)  Provides  derivative  function  f  for  ODE  y'=f(x,y). 
c  Subroutine  must  be  defined  with  parameters  F(X,Y,YP,NEqn), 

c  where  YP  is  the  derivative  at  (X,Y)  and  NEqn  is  as  defined 

c  below.  F  must  be  declared  external  in  the  calling  program, 

c  NEqn:  (integer)  Number  of  equations  in  ODE  system,  input 
c  XF:  (real*8)  Intended  termination  point  for  integration  process,  input 
c  H:  (real*8)  Step  size  used.  Output  only,  not  to  be  defined  or  changed  by 
c  the  user 

c  ATol:  (real*8)  Absolute  tolerance  for  local  truncation  error,  used  in 
c  controlling  integration  according  to  formula  described  below, 

c  input 

c  RTol:  (real*8)  Relative  tolerance  for  local  truncation  error,  used  in 
c  controlling  integration  according  to  formula 

c  ErrEst<ATol+RTol* ,  input 

c  Starting:  (logical)  Initially  set  to  .true.,  thereafter  not  to  be  changed 
c  by  the  user 

c  XOut::  (real*8)  End  point  of  last  integration  step  performed, output  only, 
c  not  to  be  changed  by  the  user 

c  NCalls:  (integer)  Number  of  function  calls  since  beginning  of  integration, 
c  Output  only,  not  to  be  changed  by  the  user, 

c  NMisses:  (integer)  Number  of  tolerance  misses  since  beginning  of 
c  integration.  Output  only,  not  to  be  changed  by  the  user. 
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c  Diag:  (logical)  Set  to  true  if  only  one  integration  step  at  a  time  is 
c  desired,  which  is  the  preferred  mode  for  study  of  integrator 

c  behavior,  and  false  if  it  is  desired  to  integrate  until 

c  interpolated  output  at  X  may  be  obtained 

c  ErrEst:  (real*8)  Maximum  norm  of  estimated  error  in  last  integration  step 
c  taken,  output  only 

c  Work:  (real*8)  Work  array  which  must  be  dimensioned  at  least  2+23*NEqn 
c  for  2nd  order  and  2+52*NEqn  for  5th  order, 

c 

c* **************** **** ****************************** ******************** 
**** 

implicit  none 

logical  Starting, Diag 

integer  NCalls,NMisses,NEqn 

real*8  H,XO,YO(NEqn),X,Y(NEqn),ATol,RTol,XF,XOut,Work(2+52*NEqn) 
real*8  ErrEst 
external  F 

call  driver(XO,YO,X,Y,F,NEqn,XF,H,ATol,RTol,Starting,XOut, 

1  NCalls,NMisses,Diag,ErrEst,Work(l),Work(2), 

2  Work(2+5  *NEqn),Work(2+ 1 0*NEqn),W  ork(2+ 1 5  *NEqn), 

3  Work(3+15*NEqn), 

4  Work(3+ 1 6*NEqn),  Work(3+ 1 7  *NEqn)  ,Work(3+ 1 8  *NEqn), 

5  Work(3+ 1 9*NEqn),W  ork(3+20*NEqn),W  ork(3+2 1  *NEqn), 

6  Work(3+22*NEqn),Work(3+27*NEqn),Work(3+32*NEqn), 

7  Work(3+37*NEqn),Work(3+38*NEqn),Work(3+46*NEqn» 
end 
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subroutine  driver(XO,YO,X,Y,F,NEqn,XF,H,ATol,RTol, Starting, XOut, 

1  NCalls,NMisses,Diag,ErrEst,H01d,YIter,FStage,YIterP, 

2  X01d,YP,Yl,Y2,Y4,YlP,Y2P,YStage, 

3  YTterN,YIterS,FStageN,Y6Der,RKWork,NordVec) 

**** 

c  driver  provides  the  essential  logic  to  control  the  integration  process  and 
c  calls  the  subprograms  needed  to  carry  out  the  necessary  calculations.  When 
c  Starting  is  set  to  .true.,  start  is  called  to  obtain  an  initial  Nordsieck 
c  vector,  step  size  and  external  stage  vector,  solver  then  is  called  to 
c  carry  out  the  first  integration  step  and  the  error  estimate  is  checked  to 
c  ensure  the  result  is  within  tolerance.  If  not  the  step  size  is  halved  and 
c  Nordsieck  vector  is  rescaled  to  produce  a  new  starting  external  stage  vector, 
c  This  process  is  repeated  until  successful.  Starting  is  then  set  to  .false, 
c  and  the  integration  process  proceeds  as  NewStep  calculates  a  new  step  size  to 
c  use.  If  Diag  is  .true,  driver  returns  dim5 1  which  returns  to  be  called  again 
c  by  the  main  program,  otherwise  at  each  iteration  either  XOut  is  at  least  X, 
c  in  which  case  Nordlnterp  produces  an  interpolated  solution  value  at  X  and  a 
c  return  is  made  to  the  main  program  for  another  call,  or  an  integration  step 
c  is  carried  out.  An  integration  step  begins  by  saving  copies  of  the  current 
c  and  previous  external  stage  vectors.  Then  rescale  is  called  to  rescale  the 
c  current  external  stage  vector,  a  DIMSIM  step  is  performed  using  a  call  to 
c  solver,  and  ErrEst  is  checked  to  be  sure  the  error  is  within  tolerance.  If 
c  not,  a  tolerance  miss  is  counted  and  the  step  size  is  halved.  Rescaling 
c  again  takes  place  using  the  saved  values,  and  this  process  continues  until 
c  the  step  is  within  tolerance  or  MaxTries  is  exceeded,  generating  an  error 
c  message  and  STOP.  For  normal  execution  final  XOut  is  XF. 
c 

c  Parameters: 

c  X0, Y0, X, Y,F,NEqn,XF,H,ATol,RTol, Starting, XOut, NCalls,NMisses, Diag, ErrEst: 
c  As  in  dim5 1 ,  see  Section  a  above. 

c  HOld:  (real*8)  Step  size  used  for  previous  successful  integration  step 
c  YIter:  (real*  8  array)  External  stage  vector  for  current  step,  dimensioned 
c  (NEqn,5) 

c  FStage:  (real*8  array)  Vector  of  derivative  values  computed  from  internal 
c  stage  values  for  most  recent  successful  step,  dimensioned  (NEqn,5). 
c  YIterP:  (real*8  array)  External  stage  vector  for  previous  successful  step, 
c  dimensioned  (NEqn, 5) 
c  XOld:  (real*8)  Previous  successful  step  point 
c  YP:  (real*8  array)  Derivative  at  X0,  =f(X0,Y0).  Dimensioned  NEqn. 
c  Y1,Y2,Y4,Y1P,Y2P,  Y6Der:  (real*8  arrays)  Storage  needed  for  subroutine  start, 
c  Each  dimensioned  NEqn. 

c  YIterN,YIterS,  FStageN:  (real  *8  arrays)  Storage  needed  for  temporary  storage 
c  of  current  and  previous  external  stage  vectors  and  current  vector  of 
c  derivatives  at  internal  stage  points  until  it  is  verified  that  step 

c  was  successful,  dimensioned  (NEqn, 5). 

c  YStage:  (real*8  array)  Storage  needed  for  subroutine  solver,  dimensioned 
c  NEqn. 

c  RKWork:  (real*  8  array)  Storage  needed  for  subroutine  start,  dimensioned 
c  8*NEqn. 

c  NordVec:  (real*8  array)  At  first  step,  Nordsieck  vector  at  X0.  Subsequently, 
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c  Nordsieck  vector  at  XOut.  Dimensioned  (NEqn,6) 

*** 

c 

c  Driver  controls  the  integration  process, 
c 

implicit  none 
logical  Starting,Diag 

real*8  H,XO,YO(NEqn),ErrEst,FStage(NEqn,5),YIter(NEqn,5),X,XOut 
real*8  YIterN(NEqn,5),YIterP(NEqn,5),Y(NEqn),X01d,H01d,HFac,delta 
real*8  Theta,Eps,YIterS(NEqn,5),HNew,YP(NEqn),ATol,RTol,YlP(NEqn) 
real*8  FStageN(NEqn,5),Y  1  (NEqn),  YStage(NEqn) 

real*8  NewStep,Xl,norm,Tol,Y2(NEqn),Y6Der(NEqn),NordVec(NEqn,6),XF 
real*8  C(5),A(5,5),B(5,5),vT(5),W(5,6),BT(6,5),Beta(5),Gamma(5) 
real*8  RKWork(8*NEqn),Y2P(NEqn),Y4(NEqn) 
common  /Method/C, A, B,vT,W,BT, Beta, Gamma 
save  Method 

integer  NCalls,NMisses,iii,jjj,NEqn,JEqn,NTries,MaxTries 
external  F 

parameter  (Eps=2.2D-16,MaxTries=10) 
if  (Starting)  then 
NCalls=0 
NMisses=0 

Tol=ATol+RT  ol*norm(  Y 0,NEqn) 
c 

c  Compute  starting  values 
c 

call  start(YIter, F,X0,Y0, NEqn, H, YP,Tol,Yl , YIP,  Y6Der,Y2, 

1  Y2P,Y4,NordVec,RKWork) 

H01d=H 

NCalls=NCalls+21 

delta=l 

c 

c  Iterate  until  H  yields  error  estimate  within  tolerance 
c 

1  continue 
c 

c  Compute  First  Step  with  Error  Estimate  and  Interp  Parameters 
c 

call  Solver(X0,NEqn,YIter,H,HOld,F,ErrEst,FStage,YIterP, 

1  Starting,YStage) 

NCalls=NCalls+5 
if  (ErrEst  .gt.  Tol)  then 

print  *, 'FIRST  STEP  MISSED!!!  H’,H,'X0=’,X0 
NMisses=NMisses+ 1 
H=H/2 
delta=delta/2 
do  JEqn=l,NEqn 
YIter(JEqn,l)=YO(JEqn) 
do  iii=2,5 

YIter(JEqn,iii)= Y  O(JEqn) 
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dojjj=2,6 

YIter(JEqn,iii)=YIter(JEqn,iii)+ 

1  W  (i  ii  ,j jj )  *N ord  V ec(  JEqn,jj j )  *  delta*  *  (jjj  - 1 ) 

end  do 
end  do 
end  do 
goto  1 
endif 

c  HNew=NewStep(H,Tol,2*ErrEst) 

XOut=XO+H 
H01d=H 
c  H=HNew 

Starting  =  .false, 
if  (Diag)  then 
do  JEqn=l.NEqn 
Y  (JEqn)= YIter  (JEqn,  1 ) 
end  do 
RETURN 
endif 
endif 
c 

c  Iterate  until  X  is  reached  or  exceeded, 
continue 

if  ((XOut  .ge.  X)  .and.  .not.  Diag)  then 
c 

c  Calculate  Y  using  Nordlnterp 
c 

Theta=(X-X01d)/H01d 

call  NordInterp(Y,H01d, Theta, FStage,YIterP,NordVec,NEqn) 
else 

NTries=0 

X01d=X0ut 

T  ol= AT  ol+RT  ol  *norm(  YIter(  1 , 1  ),NEqn) 
c 

c  Iterate  until  error  estimate  is  within  tolerance, 
c 

3  continue 

XOut=X01d+H 
if  (XOut  .gt.  XF)  then 
H=XF-X01d 
XOut=XF 
endif 

do  JEqn=l,NEqn 
do  iii=l,5 

YIterS(JEqn,iii)=YIterP(JEqn,iii) 
YIterN(JEqn,iii)=YIter(JEqn,iii) 
end  do 
end  do 
c 

c  Rescale  iteration  vector  for  new  step  size 
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call  rescale(YIterN,H,H01d,FStage,  YIterS  ,NEqn) 
call  Solver(X01d,NEqn,YIterN,H,H01d,F, ErrEst, 

1  FStageN, YIterS, Starting, YStage) 

NCalls=NCalls+5 
if  (ErrEst  .gt.  Tol)  then 
NTries=NTries+ 1 
if  (NTries  .gt.  MaxTries)  then 

print  *,'At  XOut=',XOut,',  MaxTries  exceeded’ 
print  *, Tol-, Tol 
STOP 
endif 

NMisses=NMisses+ 1 
H=H/2 
goto  3 
endif 

HNew=NewStep(H, Tol, ErrEst) 

H01d=H 
H=HNew 
do  JEqn=l,NEqn 
Y  (JEqn)= YIterN  ( JEqn,  1 ) 
do  iii=l,5 

YIterP(JEqn,iii)=YIterS(JEqn,iii) 
YIter(JEqn,iii)=YIterN(JEqn,iii) 
FStage(JEqn,iii)=FStageN(JEqn,iii) 
end  do 
end  do 

if  (.not.  Diag)  then 
goto  2 
endif 
endif 
end 
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SUBROUTINE  Solver(X0,NEqn,YIter,H,HOld,F,ErrEst,FStage,YIterP, 

1  Starting,  YStage) 

c 

c  Solver  carries  out  one  step  of  the  solution  for  an  explicit  fifth 
c  order  DIMSIM  ODE  solver  for  an  ODE  system.  U  is  assumed  to  be  the 
c  identity  matrix, 
c 

C**********************  **************  ****************************  ******* 

c  solver  applies  the  basic  DIMSIM  algorithm  to  compute  a  new  external  stage 
c  vector  and  an  error  estimate.  Different  coefficients  are  used  in  computing 
c  the  error  estimate  if  Starting  is  .true, 
c 

c  Parameters: 

c  X0:  (real* 8)  Starting  point  for  this  integration  step 
c  NEqn,H,H01d,F,ErrEst, Starting:  See  description  of  driver 
c  FStage:  (real*8)  Upon  output,  set  to  vector  of  derivatives  at  internal 
c  stage  points  for  current  step,  dimensioned  (NEqn,5) 

c  YIter:  (real*8)  Upon  input,  set  to  rescaled  previous  external  stage  vector, 
c  Upon  output,  set  to  current  external  stage  vector.  Dimensioned  (NEqn,5) 
c  YIterP:  (real*8)  Output  only,  set  to  input  value  of  YIter.  Dimensioned 
c  (NEqn,5) 

c  YStage:  (real*8)  Temporary  storage  needed  for  each  internal  stage  vector  in 

c  turn  during  calculation.  Dimensioned  NEqn. 

********************************************  *********** 

implicit  none 

real*  8  X0,H,ErrEst,YIter(NEqn,5),HOld,Delta,Err 
real*8  YStage(NEqn),YIterP(NEqn,5),FStage(NEqn,5),Eps 
real *8  C(5),A(5,5),B(5,5),vT(5),W(5,6),BT(6,5) 
real*8  Beta(5),Gamma(5),vTYIter,Fac,BetInit(5),GamInit(5) 
integer  iii,jjj,JEqn,NEqn 
logical  Starting 
parameter  (Eps=2.2d-16) 
external  F 

common  /Method/C, A,B ,vT,W,BT, Beta, Gamma 
save  Method 
Delta=H/H01d 
do  JEqn=  1  ,NEqn 
do  iii=l,5 

YIterP(JEqn,iii)=YIter(JEqn,iii) 
end  do 
end  do 

call  F(X0+C(  1  )*H, YIterP(  1 , 1  ),FStage(  1 , 1  ),NEqn) 
do  iii=2,5 
do  JEqn=l,NEqn 
Y  Stage(JEqn)= YIterP(JEqn,iii) 
do  jjj=l,iii-l 

Y  Stage(JEqn)=Y  Stage(  JEqn)+ 

1  H*A(iii,jjj)*FStage(JEqn,jjj) 

end  do 
end  do 
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call  F(XO+C(iii)*H,YStage,FStage(l,iii),NEqn) 
end  do 

do  JEqn=l,NEqn 
vTYIter=0 
do  iii=l,5 

vTYIter=vTYIter+vT  (iii)*YIterP(JEqn,iii) 
end  do 
do  iii=l,5 

YIter(JEqn,iii)=vTYIter 
do  jjj=l,5 

YIter(JEqn,iii)=YIter(JEqn,iii)+ 

1  H*B(iiijjj)*FStage(JEqn,jjj) 

end  do 
end  do 
end  do 
ErrEst=0 

if  (.not.  Starting)  then 
do  JEqn=l,NEqn 
Err=0 
do  iii=l,5 

Err=Err+H*Beta(iii)*FStage(JEqn,iii)+ 

1  Gamma(iii)*YIterP(JEqn,iii) 
end  do 

if  (abs(Err)  .gt.  ErrEst)  then 
ErrEst=abs(Err) 
endif 
end  do 

Fac=Delta**4/(18. 82229382702224d0+l. 07997453 1970991d2*Delta+ 

1  1 .43 3089043 34375d2*Delta*  *2+ 

2  5.5 1 3374496429978dl  *Delta**3+Delta**4) 

ErrEst=Fac*ErrEst 

else 

Betlnit(  1  )=0. 181 287554525493 6d0 

Betlnit(2)=0.3105254335256836d0 

BetInit(3)=-0. 3258295 1053151 99d0 

Betlnit(4)=0. 1 138267699 1 95702d0 

Betlnit(5)=-0.01051030018603054d0 

Gamlnit(  1 )=-0. 135 132430963246 1 dO 

Gamlnit(2)=0 

Gaxnlnit(3)=0 

Gamlnit(4)=0 

Gamlnit(5)=0. 1351 32430963246 1  dO 
do  JEqn=l,NEqn 
Err=0 
do  iii=l,5 

Err=Err+H*BetInit(iii)*FStage(JEqn,iii)+ 

1  GamInit(iii)*YIterP(JEqn,iii) 
end  do 

if  (abs(Err)  .gt.  ErrEst)  then 
ErrEst=abs(Err) 
endif 
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end  do 
endif 
end 
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subroutine  start(YIterP,F,XO,YO,NEqn,H,YP,Tol,Yl ,YlP,Y6Der, 

1  Y2,Y2P,Y4,YNord,RKWork) 

c*********************************************************************** 

*** 

c  start  first  calculates  the  derivative  at  the  initial  point  X0.  A  scaling 
c  factor  g  to  estimate  the  size  of  the  initial  Nordsieck  vector  is  computed  and 
c  this  is  used  to  determine  hO.  In  most  cases  g  will  be  the  smaller  of  the 
c  norm  of  the  function  and  the  first  derivative,  but  in  case  one  is  smaller 
c  than  1000*eps  (machine  epsilon)  the  larger  is  taken,  and  in  no  case  is  g 
c  taken  as  less  than  10-6.  A  6th  order  Runge-Kutta  method  is  used  by  calling 
c  RK6  to  compute  the  solution  Y4  at  X0+h0/4,  Y2  at  X0+h0/2  (from  Y4),  and  Y1 
c  at  XO+hO  (from  Y2).  These  are  then  used  to  calculate  a  Nordsieck  vector  for 
c  computation  of  the  initial  external  stage  vector,  plus  an  extra  derivative 
c  for  computation  of  the  appropriate  size  for  the  first  DIMSIM  step.  Both 
c  the  Nordsieck  vector  and  the  intial  external  stage  vector  are  returned  in 
c  case  tolerance  is  not  met. 

c  YIterP:  (real*8)  Initial  external  stage  vector,  output  dimensioned  (NEqn,5) 
c  F,X0,Y0,NEqn:  Seedim51  description 
c  H:  (real  *8)  Output  value  for  calculated  initial  DIMSIM  step  size 
c  YP:  (real* 8)  Derivative  of  the  solution  at  X0,  calculated  using  y'=F(X0,Y0), 
c  output  dimensioned  NEqn 

c  Tol:  (real*8)  Initial  tolerance  computed  as  described  in  section  a  from 
c  RTol,  Y0  and  ATol,  input 

c  Y1,Y1P,Y2,Y2P,Y4:  (real*8  arrays)  Values  computed  using  6th  order 
c  Runge-Kutta  method  for  solution  at  hO,  hO/2,  and  hO/4,  and  associated 
c  derivatives  computed  using  F,  workspace  each  dimensioned  NEqn. 
c  Y6Der:  (real*8  array)  6th  derivative  of  solution  vector  at  X0,  multiplied 
c  times  ,  workspace  dimensioned  NEqn. 

c  YNord:  (real*8)  Nordsieck  vector  at  X0  for  step  size  H,  output  dimensioned 
c  (NEqn, 6) 

cRKWork:  See  driver  description 

** 

implicit  none 

real *8  YIterP(NEqn,5),X0,Y0(NEqn),H,YP(NEqn),Tol,YPP,norm,EPS 
real*8  Y1  (NEqn),YlP(NEqn),Y2(NEqn),g,hO, delta 
real*8  C(5),A(5,5),B(5,5),vT(5),W(5,6),BT(6,5),Beta(5),Gamma(5) 
real*8  Y6Der(NEqn),YNord(NEqn,6),Y2P(NEqn),Y4(NEqn) 
real *8  RKWork(8*NEqn) 
integer  NEqn,JEqn,iii,jjj 
parameter  (EPS=2.2d-16) 
external  F 

common  /method/C, A, B  ,vT,W,BT,Beta, Gamma 
save  Method 
call  f(X0,  Y 0,  YP,NEqn) 
c 

c  Compute  a  value  for  hO 
c 

g=dmin  1  (norm(YO,NEqn),norm(YP,NEqn)) 

if  (g  .It.  1000*eps)  g=dmaxl(norm(yO, NEqn), norm(YP, NEqn)) 

g=max(g,1.0d-6) 


126 


NAWCWPNS  TP  8340 


h0=  10**(5.0d0/6)*(eps/g)**(  1 ,0d0/6) 

Compute  a  value  for  y4  at  X0+H0/4  using  sixth  order 
Runge-Kutta  method 

call  RK6(Y4,X0,Y0,h0/4,F,NEqn,RKWork(  1  +7*NEqn),RKWork(  1 ), 

1  RKWork(  1+NEqn), 

1  RKWork(  l+2*NEqn),RKWork(  1  +3  *NEqn),RKW  ork(  1  +4*NEqn), 

2  RKWork(  l+5*NEqn),RKWork(  1  +6*NEqn)) 

Compute  a  value  for  y2  and  y2P  at  X0+H0/2  using  sixth  order 
Runge-Kutta  method 

call  RK6(Y2,X0+h0/4,Y4,h0/4,F,NEqn,RKWork(  1  +7*NEqn),RKWork(  1 ), 
1  RKWork(  1+NEqn), 

1  RKW  ork(  1  +2*NEqn),RKW  ork(  1+3  *NEqn)  ,RKW  ork(  1  +4*NEqn), 

2  RKWork(  1  +5*NEqn),RKWork(  1  +6*NEqn)) 
call  F(X0+h0/2,Y2,Y2P,NEqn) 

Compute  a  value  for  yl  and  ylP  at  X0+H0  using  sixth  order 
Runge-Kutta  method 

call  RK6(Y1  ,X0+h0/2,Y2,h0/2,F,NEqn,RKWork(  1  +7*NEqn),RKWork(  1 ), 
1  RKW  ork(  1  +NEqn) , 

1  RKWork(  1  +2*NEqn),RKWork(  1  +3*NEqn),RKWork(  1  +4*NEqn), 

2  RKW  ork(  1+5  *NEqn)  ,RKW  ork(  1  +6  *NEqn) ) 
call  F(X0+h0,  Y 1  ,Y  1  P,NEqn) 

Compute  sixth  derivative  times  hOA6 

do  JEqn=l,NEqn 

Y6Der(JEqn)=1280*(-90*Y0(JEqn)-22*Yl(JEqn)-144*Y2(JEqn)+ 

1  256*Y4(JEqn)+H0*(-9*YP(JEqn)+ 

2  3*Y  lP(JEqn)+36*Y2P(JEqn))) 
end  do 

Compute  initial  step  size 

H=hO*(Tol/(4*(5.5 1 8667640362375d-5)*norm(Y6Der,NEqn))) 

1  **(1.0d0/6) 

if  (H  .It.  100*EPS)  then 

print  *,'START-Tolerance  too  tight  or  H  too  small' 
print  *,'h0=',h0,'H=',H,'New  H=1.0d-8' 

H=1.0d-8 

endif 

Compute  starting  vector 

delta=H/hO 
do  JEqn=l,NEqn 
YNord(  JEqn,  1  )= YO(JEqn) 
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YNord(JEqn,2)=H*YP(JEqn) 

YNord(JEqn,3)=2*(-567*Y0(JEqn)-25*Yl(JEqn)-432*Y2(JEqn)+ 

1  1 024*  Y4(JEqn)+HO*(-90*  YP(JEqn)+3  *  Y 1  P(  JEqn)+ 

2  72*Y2P(JEqn)))/9 
YNord(JEqn,3)=YNord(JEqn,3)*delta**2 
YNord(JEqn,4)=1836*Y0(JEqn)+148*Yl(JEqn)+21 12*Y2(JEqn)- 

1  4096*Y4(JEqn)+H0*(222*YP(JEqn)- 

2  1 8  *  Y 1  P(  JEqn)-3 84*  Y 2P(  JEqn)) 
YNord(JEqn,4)=YNord(JEqn,4)*delta**3 

YNord(JEqn,5)=32*(-1323*Y0(JEqn)-169*Yl(JEqn)-1836*Y2(JEqn)+ 

1  3328*  Y4(JEqn)+H0*(- 144*  YP(JEqn)+ 

2  21  *YlP(JEqn)+378*Y2P(JEqn)))/3 
YNord(JEqn,5)=YNord(JEqn,5)*delta**4 
YNord(JEqn,6)=160*(378*Y0(JEqn)+70*Yl(JEqn)+576*Y2(JEqn)- 

1  1024*Y4(JEqn)+H0*(39*YP(JEqn)- 

2  9*YlP(JEqn)-132*Y2P(JEqn))) 
YNord(JEqn,6)=YNord(JEqn,6)*delta**5 

end  do 


c 

c  We  make  use  of  the  special  form  of  W  for  c0=0,  explicit  methods 
c 

do  JEqn=l,NEqn 
YIterP(JEqn,  1  )= YO(JEqn) 
do  iii=2,5 
YIterP(JEqn,iii)=0 
do  jjj=l,6 

YIterP(JEqn,iii)=YIterp(JEqn,iii)+ 

1  W  (iiijjj)  *  YNord(JEqnjjj) 

end  do 
end  do 
end  do 
end 
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subroutine  rescale(YIter,H,H01d,FStage,YIterP,NEqn) 

He********************************************************************** 

rescale  takes  a  previous  external  stage  vector  computed  using  a  step  size 
HOld  and  uses  the  DIMSIM  rescaling  algorithm  to  rescale  it  to  become  the 
appropriate  previous  external  stage  vector  for  continuing  the  computation 
with  step  size  H. 

Parameters: 

YIter:  (real*8)  Output  rescaled  extermal  stage  vector,  dimensioned  (NEqn,5) 

H:  (real*8)  Input  new  step  length 

HOld:  (real* 8)  Input  step  length  for  previous  successful  step. 

FStage:  (real*8)  Input  vector  of  derivatives  at  stage  points  for  last  step 
completed 

YIterP:  (real*8)  Input  previous  external  stage  vector,  calculated  using  step 
length  HOld,  dimensioned  (NEqn,5) 

NEqn:  Seedim51 

He********************************************************************** 

implicit  none 

real*8  Yiter(NEqn,5),H,H01d,FStage(NEqn,5),YIterP(NEqn,5),Delta 
real *8  C(5),A(5,5),B(5,5),vT(5),W(5,6),BT(6,5),Suml  ,Sum2 
real*8  Beta(5),Gamma(5) 
integer  i,j,k,JEqn,NEqn 

common  /Method/C, A,B ,vT,W,BT, Beta, Gamma 
save  Method 
Delta=H/H01d 
do  JEqn=l,NEqn 
do  i=l,5 
Sum  1=0 
Sum2=0 
do  j=l,5 

Sum2=Sum2+vT  (j )  *  YIterP( JEqn  ,j ) 
end  do 
doj=l,5 
do  k=l,6 

Suml=Suml+W(i,k)*Delta**(k-l)*BT(k,j)*FStage(JEqn,j) 
end  do 
end  do 

YIter(JEqn,i)=H01d*Suml+Sum2 
end  do 
end  do 
end 
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real*8  function  NewStep(H,Tol, ErrEst) 

0  ^  jf*  jjc  5jv  jfc  ?fC  5$C  jji  ^fc  jjc  5jC  3|C  ijc  !jc  ^  5jC  5jC  5j>  5jC  #{c  5jc  #jc  #Jc  5fc  5{C  5jC  ijC  j|l  ^  ^  5fc  5jc  ^  2»Jc  5^C  j|c  5jc  5jC  #fC  ij*  5|C  jjC  SfC  5jC  jjC  ?jc  5|C  5[C  2jc  jjc  *|C  5jc  >|C  5{C 

** 

c  The  error  estimate  is  compared  with  the  tolerance  to  produce  a  step  change 
c  factor  modifying  the  step  size  in  preparation  for  the  next  step.  The  local 
c  truncation  error  is  assumed  to  be  proportional  to  the  6th  power  of  the  step 
c  size.  The  largest  possible  step  change  is  a  factor  of  2.  A  step  change 
c  safety  factor  equal  to  .75  was  used  for  DIMEX5  so  that  the  step  size  used 
c  was  .75  of  optimal.  A  larger  factor  of  .9  was  used  successfully  for  DIMEX2. 
c 

c  Parameters: 

c  H:  (real*8)  At  input  contains  old  step  size.  At  output  contains  new  step 
c  size 

c  Tol:  (real*8)  Input  tolerance  derived  from  RTol,  Y  and  ATol  as  described 
c  in  section  a  above 

c  ErrEst:  (real*  8)  Input  current  local  error  estimate 

*1*  vL*  *2*  ^  *1*  -a*  '*!>*’  4*  «sL»  *1*  4^  *4"  *1*  *1*  4*  4>  4*  4*  4*  4*  4*  *4*  *4*  *4*  4*  44*  *4*  4*  *4^  4*  *4*  4*  4*  4*  4f  4*  *4*  4/  4f  4^  4*  4*  4?  4*  4f  4*  ^  4*  4  4?  *4f  4f  4  4^  *4*  "4* 

^  ^  ^  4*  'T'  *T*  “T*  *T*  *T*  *T*  *T*  “T*  *T*  *r*  *T*  <iv  *T*  *T*  4'  *T*  *1*  'V  'T*  *1*  'T*  *1*  'T*  *V*  "T*  ■’T*  *1*  'I'  *T*  'T*  *1*  'r  'P  ■'T*  *T*  *T“ 

** 

implicit  none 

real*8  ErrEst, RTol,ATol,Tol,HFac,H,Eps,YNorm 
parameter  (Eps=2.2d-16) 

HFac=.75d0*(Tol/abs(ErrEst))**(1.0d0/6) 

HFac=dmin  1  (HFac,2.0d0) 

NewStep=HFac*H 

end 
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subroutine  NordInterp(Y,H,Theta,FStage,YIter,NordVec,NEqn) 

*  ****:{::*:  *****:*:*;{:*************************************************  ******* 


The  Nordsieck  vector  NordVec  is  calculated  at  the  current  value  of  XOut, 
using  the  previously  calculated  external  stage  vector  YIter,  the  current 
step  size  H,  and  the  current  vector  of  derivatives  at  stage  points,  FStage. 
Theta  then  is  used  to  rescale  the  Nordsieck  vector  to  correspond  to  the 
desired  output  point  and  a  5th  order  Taylor  series  approximation  Y  is  formed. 

Parameters: 

Y:  (real*8  vector)  Output  approximation  to  the  solution  at  the  desired 
output  point,  length  NEqn. 

H:  (real*8)  Input  step  size  for  last  successful  integration 
Theta:  (real*  8)  Input  between  0  and  1,  ratio  (X-X01d)/H,  where  X  is  the 
desired  output  point 

FStage:  (real*8  vector)  Input  vector  of  derivatives  at  stage  points  for 
last  successful  integration,  dimensioned  (NEqn, 5) 

YIter:  (real*8  vector)  Input  previous  external  stage  vector,  dimensioned 
(NEqn,5) 

NordVec:  (real*8  vector)  Workspace  (output  for  waveform  relaxation)  for 
Nordsieck  vector  at  current  XOut  and  last  step  size  used. 

NEqn:  Seedim51 


implicit  none 

real*8  H,Theta,FStage(NEqn,5),YIter(NEqn,5),Y(NEqn) 
real*8  NordVec(NEqn,6),Bt(6,5),vT(5),vTIter,del,C(5),A(5,5) 
real*8  B(5,5),W(5,6),Beta(5),Gamma(5),Factorial 
integer  JEqn,NEqn,iii,jjj 
common  /Method/C, A,B,vT,W,Bt, Beta, Gamma 
save  Method 
do  JEqn=l,NEqn 
vTIter=0 
do  iii=l,5 

vTIter=vTIter+vT(iii)*YIter(JEqn,iii) 
end  do 
do  iii=l,6 

NordVec(JEqn,iii)=vTIter 
do  jjj=l,5 

NordVec(JEqn,iii)=NordVec(JEqn,iii)+ 

1  H*Bt(iii,jjj)*FStage(JEqn,jjj) 

end  do 
end  do 

del=-(l  -theta) 

Y  (JEqn)=NordVec(  JEqn,  1 ) 

factorials 

do  iii=2,6 

Y ( JEqn)= Y ( JEqn)+NordV ec(JEqn,iii)  *del  *  *  (iii- 1  )/Factorial 
Factorial=Factorial*iii 
end  do 
end  do 
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132 


NAWCWPNS  TP  8340 


subroutine  RK6(Y,X0,Y0,H,F,NEqn,YP,YP2,YP3,YP4,YP5,YP6,YP7,YTemp) 

£**:*::*::*:*:*:*;$:*  ***:}:  ***************************************************  ****** 

*** 

c  A  simple  6th  order  Runge-Kutta  solver  with  7  stages,  used  to  provide  solution 
c  valuesneeded  for  starting  method.  Solution  Y  for  initial  values  (X0,Y0)  is 
c  found  at  XO+H  for  problem  Y'=F(X,Y).  The  method  is  described  in  Butcher, 
c  pp.  203-205. 
c 

c  Parameters: 

c  Y:  (real*8  vector)  Solution  at  XO+H,  dimensioned  NEqn. 
c  X0:  (real*8)  Initial  value  of  independent  variable. 

c  Y0:  (real*8  vector)  Initial  value  of  dependent  variable,  dimensioned  NEqn 
c  H:  (real*8)  step  size  used 
cF:  seedim51 
cNEqn:  seedim51 

c  YP,YP2,YP3,YP4,YP5,YP6,YP7:  (real* 8  vectors)  Workspace  used  for  storing 
c  derivatives  of  stages,  each  dimensioned  NEqn. 

c  YTemp:  (real*8  vector)  Workspace  used  for  storing  stage  vectors 
c  successfully,  dimensioned  NEqn 

^^^^^** ************************************************* **************** 
*** 

implicit  none 
integer  JEqn,NEqn 

real *8  H,X0,Y0(NEqn),YP(NEqn),YP2(NEqn),YP3(NEqn),YP4(NEqn) 
real* 8  Y (NEqn) , YP5 (NEqn), YP6(NEqn) , YP7 (NEqn) , YT emp(NEqn) 
external  F 

call  F(X0,Y0,YP,NEqn) 
do  JEqn=l,NEqn 

YTemp(JEqn)=Y0(JEqn)+H*YP(JEqn)/3 
end  do 

call  F(X0+H/3, YTemp,  YP2, NEqn) 
do  JEqn=l,NEqn 

YTemp(JEqn)=Y0(JEqn)+2*H*YP2(JEqn)/3 
end  do 

call  F(X0+2*H/3, YTemp, YP3, NEqn) 
do  JEqn=l,NEqn 

YTemp(JEqn)=Y0(JEqn)+H*(YP(JEqn)+4*YP2(JEqn)-YP3(JEqn))/12 
end  do 

call  F  (XO+H/3 ,  YTemp,  YP4,NEqn) 
do  JEqn=l,NEqn 

YTemp(JEqn)=Y0(JEqn)+H*(25*YP(JEqn)-110*YP2(JEqn)+ 

1  35*YP3(JEqn)+90*YP4(JEqn))/48 

end  do 

call  F(X0+5*H/6, YTemp, YP5, NEqn) 
do  JEqn=l,NEqn 

YTemp(JEqn)=Y0(JEqn)+H*(18*YP(JEqn)-55*YP2(JEqn)- 
1  1 5  *  YP3  (JEqn)+60*  YP4(  JEqn)+ 1 2*  YP5  ( JEqn))/ 1 20 

end  do 

call  F(X0+h/6,YTemp,YP6,NEqn) 
do  JEqn=l,NEqn 

YT  emp(  JEqn)= Y  0(  JEqn)+H*  ((-26 1  d0/260)*  YP(JEqn)+ 
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1  (33d0/13)*YP2(JEqn)+(43d0/156)*YP3(JEqn)+ 

2  (-1 18d0/39)*YP4(JEqn)+(32d0/195)*YP5(JEqn)+ 

3  (80d0/39)*YP6(JEqn)) 
end  do 

call  F(X0+h,YTemp,YP7,NEqn) 
do  JEqn=l,NEqn 

Y(JEqn)=Y0(JEqn)+H*(13*YP(JEqn)+55*YP3(JEqn)+55*YP4(JEqn)+ 
1  32*  YP5(JEqn)+32*  YP6(JEqn)+ 1 3  *  YP7  (JEqn))/200 

end  do 
END 
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t 


block  data  Matrices 

c  The  matrices  and  vectors  defining  the  particular  method  are  set  to 
c  appropriate  values  using  data  statements  and  stored  in  the  common  region 
c  /Method/. 

q  j|C  sjc  3jC  5jC  5jC  5jC  5j%  5$C  jfC  5j?  ^  5jC  5{C  lj%  jfj  5fC  5^C  JjC  5jC  3jC  5jC  5(5  5jC  »jc  )|C  5jC  2jC  5jC  «jc  ^  3jC  5jC  *j€  jjc  5jC  5jC  «jC  #jc  #|C  5jC  5j%  3jC  #jC  5^  ^ 


real* 8  C(5),A(5,5),B(5,5),vT(5),W(5,6),BT(6,5),Beta(5),Gamma(5) 
common  /Method/C, A, B,vT,W,BT, Beta, Gamma 
save  Method 

data  VT/-.2406956 1 553862 1 5d0, 1 .260494575847145  ldO, 

1  -2.48 1 2693924523267 dO,  1 .9 1 99070083032958d0, 

2  5.4 1 5634238405067d- 1/ 
data  C  /0,2.5d-l,5.0d-l,7.5d-l,l/ 

data  A/0, 1 . 176528 1703 106688d0, 1 .9805793463233 1 9 1  dO, 

1  3 .0532835 108392395d0, 1 .4 1 93325269467698d0,2*0, 

2  .4017181 027378085d0,  .734996 1462028882d0, 

3  2.6534897473 12533  Id0,3*0,.2672357626475791d0, 

4  -2.2778532945468265d0, 4*0,1. 1978905088172778d0, 5*0/ 
data  B/3.163023914364555736d0,3.250176692142333514d0, 

1  3.508528033848969459d0,4.176223954923772036d0, 

2  3.008220785304765914d0, 1.9743392460505406810, 

3  1 .53 1 978 1 3493942957d0,0.327374204 1 84027625d0, 

4  -2.6 1 827 1 7 1 9393 1 1 992d0,2. 11453341 958770507d0, 

5  -0.81 0425 120055628 1 3d0, 0.0979082 1 3277705203d0, 

6  2.2390605 19232953536d0,7.03900409877443037d0, 

7  0.3572048837825639d0,0.540922018802018401d0, 

8  -0.422272425642426043d0, -2.097452509375452 1 55d0, 

8  -5.2884360132659356d0, 

9  -1.90425330857598604d0, 0.0550786541963 1683844d0, 

1  -0.46 138007 1 6699075 17  Id0,-0.936868983593822539d0, 

2  - 1 .69 1 09702737 1 050 1 62d0,-0.64562655527099952d0/ 
data  Bt/3 . 1 630239 1 4364555d0,0, 1 , 1 ,46666666666667d  1 ,96,256, 

1  1.974339246050541d0,0,-5.333333333333333d0, 

2  -7 .466666666666667 d  1  ,-448,- 1 024,-0. 8 1 0425 1 20055628d0,0, 

3  1 2, 1 52,768, 1 536,0.5409220 1 88020 1 84d0,0,- 1 6, 

4  - 1 .386666666666667d2,-576,- 1 024,0.055078654 1 963 1 683dO, 

5  l,8.333333333333333d0,4.666666666666667dl, 160,256/ 
data  W/l ,  1 , 1 , 1 , 1 ,0,-0.926528 1 703 1 0669d0,- 1 .88229744906 1 1 28d0, 

1  -3.3055 15419689707 dO, -1.992859488529754d0,0, 

2  0.03 1 25 dO, 0.02457 04743 15547 88d0,-0.03 611 69 178745 11 6d0, 

3  0.077 1 3 632883232 1 62d0,0,0.002604 1 66666666666d0, 

4  0.00827 964262277682d0,0.0 1 393 9400 1 002 1 236dO, 

5  0.03 1 57006827 6643 94d0, 0,0.000 1 627604 1 66666666d0, 

6  0.00 1558025774 1 2029 1 dO, 0.005702 1 29564 1 054 1 5dO, 

7  -0.0020148623151 15681d0,0,8.13802083333333d-6, 

8  0.0001 950328608825 182d0, 0.001 161984318267553d0, 

9  -0.001959141967572072d0/ 

data  Beta/-2.385388840413769d3,-1.706352417381188d2, 

1  1.870907224366729d2,-2.323292678834744dl,0/ 
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data  Gamma  /1 .268306923377 1 63d3,0,- 1 .22457220 1 88 1436d3,0, 
1  -4.373472 149572754dl/ 

end 
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